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1. Introduction 


The mathematical formulation of the engineering optimization problem is 

min XU)) 

subject to £i({x))^0, i=l,q 


( 1 ) 


where 

{x} is an nxl matrbc of design variables, 
f({x}) is the objective function, and 
gj({x}) are constraint equations. 

Evaluation of the objective function and constraint equations in Equation (1) can be very 
expensive in a computational sense. Thus, it is desirable to use as few evaluations as 
possible in obtaining its solution. In solving Equation (1), one approach is to develop 
approximations to the objective function and/or restraint equations and then to solve 
Equation (1) using these approximations in place of the original functions. These 
approximations are referred to as response surfaces. 

The desirability of using response surfaces depends upon the number of functional 
evaluations required to build the response surfaces compared to the number required in the 
direct solution of Equation (1) without approximations. The present study is concerned with 
evaluating the performance of response surfaces so that a decision can be made as to their 
effectiveness in optimization applications. In particular, this study focuses on how the 
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quality of approximations is effected by design selection. Polynomial approximations and 
neural net approximations are considered. 


To provide the groundwork for future discussion, this introductory section discusses: 

1. measures of quality of fit at the designs and measures of quality of fit over a region of 
interest and 

2. the methodology used to build the approximations. 

1.1 Quality of Fit 

Let us consider a problem with n design variables, the components of the vector {x} = {x l5 
x 2 ,.. JCn}*. A total of N designs will be considered: {x}j, j = 1,N. At the designs {x}j, let 
y. = the value of the function to be approximated and 
yj = the value of the approximating function. 

The approximating function, y, should closely match the function, y, not only at the designs, 
{x}j, but over the entire region of interest. 

1.1.1 Fit at the designs 

The approximating function y closely approximates the function y when s is small where 



and where S 2 is the sum of the squares of the residuals thus 


2 



( 3 ) 


1 

Let y be the average value of the designs, Thus 

. f»< ( 4) 



In this study, one measure of the closeness of fit to be considered is the non-dimensional 
value v where 


v.1 


ZXyrfi) 2 

i 

N 


* 100 


(5) 


The coefficient v is the non-dimensional root mean square (RMS) error at the designs. 
Thus, v = 0 is a necessary and sufficient condition that the approximating function fit the 
actual function at the N design points. 


1.1.2 Overall fit 

Just because the approximating function exactly fits the function at N designs does not 
guarantee that it gives a good fit over the region of interest. It is therefore desirable over 
the region of interest to have a measure of the quality of overall fit. Several examples of 
this study considers a two dimensional region of interest. For these problems, the 
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rectangular region of interest is overlaid with a 31x31 evenly spaced grid of points. The 
value of the function and the approximating function is then compared at these NG = 961 
evenly spaced grid of points. Other examples consider a rectangular n dimensional region 
of interest. These regions of interest are also overlaid with a evenly spaced grid of points. 
The value of the function and the approximating function are then compared at these NG 
grid points. For these examples, a measure of the quality of overall fit is taken as 


V G= 


1 


NG 


Z(y r yfiNG 
-1 * 100 




( 6 ) 


where y G is the average value of y at the grid points. A small value of v G indicates that the 
approximating function did a good job of approximation over the region of interest. 

1 2. Polynomial Approximations 

With the polynomial response surface approach, the approximating function is taken as an 
m=k+l term polynomial expression [1-3] thus 

y=b 0 +b l X l +...b t X k (?) 

where Xj is some expression involving the design variables. For example, a second order 
polynomial approximation in two variables could be of the form 
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y=K +b l X l +b 2 X 2 +b 3 X l + K X l X 2 +b 5 X i (8) 

The value of the function to be approximated at the N designs can be used to determine the 
m=k+l undetermined coefficients in the polynomial expression. For the N designs, 
Equation (7) yields 


Vi 



K 



1 X h ... 

b l 





y N . 


1 X. ... x. 

l N 

\ 


or 


[Y)=[Z\{ b ) 


( 10 ) 


where {Y} is an Nxl matrix, [Z] is an Nxm matrix, and {b} is an mxl matrix. 

1.2.1 Exactly-determined approximation 

When N = m, the approximation is exactly-determined and the matrix {b} can be determined 
from Equation (10). 

1 .2.2 Over-determined approximation 

With N>m, Equation (10) can be solved in a least squares sense thus [1-3] 
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\zf\Y\-izmw 


(ii) 


or 


{b}=amz\r 1 [z\ t {Y} 


( 12 ) 


Equation (12) in effect, chooses the terms of {b} so as to minimize the square of the 
residual as defined in Equation (2). 

1.2.3 Under-determined approximation 

When N<m, the approximation is under-determined. A solution can be obtained by 
choosing the terms of {b} so as to minimize the square of the residual as defined in 
Equation (2). However, a direct solution can be obtained by using the concept of pseudo- 
inverse [4,5]. Assume that the rank of matrix [Z] is N and define the pseudo-inverse of 
matrix Z, Z* thus 

[ZT=[ZT([Z][Z]r l (13) 


where t denotes transpose. Solution of Equation (10) is then 

[b}=[Z]'{Y}+[Q]{w) 


(14) 
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where {w} is an (m-N) column matrix of arbitrary coefficients and [Q] is a mx(m-N) matrix 
formed from any m-N independent columns of the matrix [R] thus 




(15) 


One solution to Equation (14) is to take all the arbitrary terms of {w} as zero giving 

{b)=[Z]'{Y} (16) 

The basic solution to Equation (10) is Equation (16). Using that equation, at the designs, 
{x}j, the value of ^ matches the value of yj. If w ; is the ith term in matrix {w} and {q>, is 
the ith column of matrix [Q], then at the designs, {x}j, ^=0 when 

(17) 


Thus, the last term of the right hand side of Equation (14) gives ^ values which match y^ 
at the designs, {x}j, for any values of Wj. 

1.3 Artificial Neural Nets 

While the initial motivation for developing artificial neural nets was to develop computer 
models that could imitate certain brain functions, neural nets can be thought of as another 
way of developing a response surface. Different types of neural nets are available [6,7], but 
the type of neural nets considered in this paper are back propagation nets with one hidden 
layer as shown in Figure 1. This type of neural net has been used previously to develop 
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response surfaces [8-12] and is capable, with enough nodes on the hidden layer, of 
approximating any continuous function [13]. 


For the neural net of Figure 1, associated with each node on the hidden layer, node j, and 
each output node, node k, are coefficients or weights, 0j and 0 k , respectively. These weights 
are referred to as the biases. Associated with each path, from an input node i to node j on 
the hidden layer, is an associated weight, Wy and from node j on the hidden layer to output 
node k is an associated weight Wj k . Let ^ be inputs entered at node i. Node j on the 
hidden layer receives weighted inputs, Wyq,. It sums these inputs and uses an activation 
function to yield an output Tj. The activation function considered in this paper is the 
sigmoid function [6,7] 


r _ 1 


( 18 ) 


Output node k then receives inputs w jk rj which are summed and used with an activation 
function to yield an output s k . Some variation of the delta-error back propagation algorithm 
[6,7] is then used to adjust the weights on each learning try so as to reduce the values 
between the predicted and desired outputs. In this investigation, studies were performed 
using the program NEWNET [14] which was developed especially for this investigation. 
NEWNET minimizes the sum of the squares of the residuals in Equation (2) with respect 
to the weights and biases of the net. Training of the net is thus formulated as an 
unconstrained minimization problem. Solution of this minimization problem is performed 
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using the method of Davidon, Fletcher, and Powell [15-16], That algorithm performs a 
series of one dimensional searches along search directions. Search directions are 
determined by building an approximation to the inverse Hessian matrix using gradient 
information. Gradients required by that algorithm are obtained using back-propagation. 
One-dimensional searches are performed along the search directions using an interval 
shortening routine. 
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2. Levels of Designs 


2.1 Tavlor Series Approximation 

The overriding factor which affects the accuracy of an approximation is the levels of the 
design parameters considered. It is instructive to consider a problem in two design 
variables. Suppose we wish to make a quadratic approximation of a function thus: 

y= b o +b i x i +i 2*2 + Vi + Vi*2 + V2 • • • (19) 

Consider that the exact function is evaluated at 6 design points and the information thus 
generated will be used to determine the 6 undetermined coefficients in Equation (19). 
Design variables at these design points are taken from the following sets: 

x, from the set [x. x, ..jc, } 

lJ 1 h h y (20) 

x 2 from the set [x 2j 


Here p discrete values are considered for Xj and q discrete values are considered for x 2 . 
The variable x x is said to have p levels and x 2 is said to have q levels. The problem is to 
determine the minimum levels of the design variables, p and q, required to build the 
quadratic approximation. In this regard, it is instructive to consider a Taylor series 
approximation [17] of the function about the point {Xj=0, x 2 =0}: 

y=y(0,0)+ { vy(0,0)}‘{Ax} + {ax}'[//(0,0)]{ax} + ... ( 21 ) 
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where 


{^}=[(* r o) (* 2 - 0 ) 1 ^ * 2 r 


( 22 ) 


{ yy(0,0) } = [( 


fo(o,o) fo(o,o) 

dx x dx 2 


V 


(23) 


[^(0,0)]= 


^0,0) ay 2 (0,0) 

dxf &1&2 

ay 2 (Q,Q) ^y(Q,Q) 

dx { dx 2 dx 2 


(24) 


Entering Equations (22), (23), and (24) into Equation (21) gives 


y=y( 0,0)+ 


cW0,0)_ . 3y(0,0). _ #K0,0) 2. 


dx. 


- x i + 


dx~ 


dx. 


^(0,0) ^0,0) 2 

z *i* 2 — *2 


dx x x 2 


(25) 




The derivatives in Equation (25) can be determined by finite difference equations [18]. The 
second derivative of y with respect to x } can be obtained using information at points 
indicated in Figure 2 by solid circles, the second derivative of y with respect to x 2 can be 
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obtained using information at points indicated by unfilled circles, and the mixed derivative 
can be obtained using information at points indicated by unfilled squares. 

It can be seen in Figure 2 that at least three levels of both x x and x 2 must be used to obtain 
a quadratic approximation. If three levels are not provided, not information is available to 
calculate the higher derivatives in Equation (25). A complete 3 factorial design does not 
have to be used-only 6 selected points from the complete 3 factorial design. Information 
at those 6 points allow the undetermined coefficients to be exactly determined. 

Consider now the design of Figure 3 which are also taken from the 3 factorial design. Even 
though 6 design points are used, this set of design points does not allow an approximation 
containing the x 2 2 term of Equation (25). However, with the design of Figure 3, an 
approximation of the form of Equation (26) could be obtained thus: 

y = VVi + V2 + Vi 2+ Vi*2 (26) 

With the design of Figure 3, if a solution is attempted using Equations (19) and (12), a 
singular coefficient matrix will be encountered. A solution could be attempted using the 
pseudo-inverse concept of Equations (13) and (14). However, recent studies [19] have 
shown that non-unique solutions are obtained with this technique. Non-uniqueness makes 
these solutions undesirable. Using Equations (26) and (12), a slightly over-determined 
approximation is obtained. 
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Recent studies have found that the numerical performance of neural network 
approximations and polynomial approximations with the same number of associated 
undetermined parameters is comparable [19]. Thus, it is not expected that neural nets as 
approximators will perform better than polynomials when there are inadequacies in the 
training design, as in Figure 3. The next example investigates performance of both 
polynomial and neural net approximations. 


2.2 Example 

Consider the function 


y = 1 + jc , +x 2 +x 3 +xi+x 1 x z +x ] x 3 +x% +x^c 3 +x% 


(27) 


In the first phase of the investigation, approximations are to be made of this function using 
the design of Figure 4. The star pattern of design points in Figure 4 does not allow mixed 
derivatives of the function to be calculated using finite difference type formulae but does 
permit the other second derivatives to be calculated. Thus, information is available to make 
a polynomial approximation of the form 

y = W i + V2 + *wVi 2 + Vi + V? (28) 


The function y was evaluated at the design points shown in Figure 4 yielding 7 training pairs 
for calculating the 7 undetermined parameters in Equation (28). The value of the 
approximating function y was then evaluated at a 5x5x5 grid of designs. These values of y 
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were then used to evaluate v G from Equation (6). The value of v G obtained is shown in the 


first line of Table 2.1. 

Table 2.1. Performance of Approximations for Various Designs 


Number 

Designs 

Points 

Description 

Polynomial 

Approximation 

Neural Net 
Approximation 

No. 

Para. 

s 

> — ✓ 

o 

> 

ih 

No. 

Para. 

No. 

Apx. 

v G (%) 

7 

Star--see 
Figure 4 

7 

34.6 

2 

11 

10 

25.5-97.3 

12 

Star-see 
Figure 5 

7 

34.6 

2 

11 

10 

32.9-93.5 

10 

Computer 

Generated 

10 

0.0 

2 

11 

10 

36.6-36.9 

3 

16 

10 

21.9-36.7 

27 

3 factorial 

i 

10 

0.0 

3 

16 

2 

16.6-16.7 

n 

21 

2 

16.6-16.9 

125 

5 factorial 

10 

0.0 

8 

41 

1 

3.7 


A neural net approximation was then considered. Previous studies [19] have indicated that 
it is desirable to have more training pairs than the number of undetermined parameters 
(weights and biases) associated with the net. If fewer training pairs than undetermined 
parameters are used, non-unique approximations should be expected. For a neural net with 
one hidden layer as shown in Figure 1, there are 6 parameters associated with a net with 
one node on the hidden layer and 11 parameters associated with a net with two nodes on 
the hidden layer. It was considered that one node on the hidden layer would yield an 
inadequate approximation. Thus 2 nodes on the hidden layer were considered. Thus, the 
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neural net approximation is under-determined. That is to say that there are fewer training 
pairs than there are undetermined parameters associated with the approximation. Non- 
unique approximations are to be expected. Indeed, this was the case. The 8 training pairs 
were used to make 10 different approximations by having training commence from a 
different randomly selected set of weights and biases. Once the nets were trained, the value 
of the approximating function, y, was generated at the 5x5x5 set one grid points and the 
value of v G was developed. The range of the values obtained is shown in Table 2.1. One 
can see that a large range of values is obtained. The best neural net approximation is only 
slightly better than the polynomial approximation while the worst neural net approximation 
is considerably worse. Just as with the polynomial approximation, the designs used to train 
the approximation can not yield information necessary to capture essential features of the 
function to be approximated. 

The 12 designs of Figure 5 were next used in the training of a polynomial approximation 
and a 2 node neural net approximation. Even though more designs are used here than in 
Figure 4, the additional designs selected do not yield any more information about the nature 
of the function being approximated. Information is still not available for determining the 
mixed derivatives of the function to be approximated. Thus, the polynomial approximation 
of Equation (26) was considered. As there are now more training pairs than there are 
undetermined parameters, the approximation obtained is over-determined. As no new 
information is available with the 12 designs, the same polynomial approximation and thus 
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the same v G as before are obtained. The value of v G is shown in the second line of Table 

2 . 1 . 

A neural net with 2 nodes on the hidden layer was then trained with the 12 training pairs. 
The net was trained 10 times starting from different randomly selected sets of weights and 
biases. Even thought the number of training pairs, 12, is greater than the number of 
undetermined parameters associated with the net, 11, non-unique approximations were 
obtained as can be seen in Table 2.1. Thus, it can be concluded that for neural net 
approximations, having more training pairs than the number of asso ciated undetermined 
parameters is only a necessary condition for obtaining a unique approxim ation but that it 
is not a sufficient condition . As the 12 designs offered no new information about the 
function being approximated over that offered by the 8 designs, then just as with the 8 
design case, non-unique approximations were obtained. 

The program DESIGNS [20], which was developed for this project, was used to generate 10 
designs which contain the information necessary for calculating the 10 undetermined 
coefficients of the complete quadratic approximation of the form: 

y=b 0 +b l x l +b 2 x 2 +b 3 x 3 +b 4 x 1 +b s x 2 +b 6 x 3 +b 7 x l x 2 +b t p 1 x 3 +b^x 2 x 3 ^9) 

The location of these design points is shown in Figure 6. The polynomial approximation 
obtained by training the polynomial of Equation (29) with the computer generated designs 
exactly duplicated the test function of Equation (27). Thus, v G for the 5x5x5 grid of points 


16 



was zero as seen in the third line of Table 2.1. 


A neural net with 2 nodes on the hidden layer with 6 associated undetermined parameters 
and a neural net with 3 nodes on the hidden layer and 11 associated undetermined 
parameters were then trained 10 times with the computer generated training pairs. Each 
training started from a different randomly selected set of weights and biases. For the case 
of 2 nodes on the hidden layer, the approximation generated was over-determined and a 
unique approximation was obtained (the small range of v G obtained most likely results from 
the exit criteria employed in the training algorithm). For the case of 3 nodes on the hidden 
layer, there are 11 associated undetermined parameters but only 10 training pairs. Thus the 
approximation is under-determined and a non unique approximation is obtained as can be 
seen in Table 2.1. 

The performance of the neural net approximations was much poorer than that of the 
polynomial approximation on this problem. This poorer performance may be in part 
because the problem is biased towards the polynomial approximation as the function being 
approximated is 2 second order polynomial. 

A complete 3 3 factorial design and a 5 3 factorial design were considered to see if good 
results could be obtained with the neural nets if more training pairs were employed. Indeed 
this was the case. However, many more training pairs were required to get a good 
approximation than were required with the polynomial approximation. The extra training 
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pairs were wasted on the polynomial approximation. Ten correctly selected training pairs 
is all that is required to get an exact second order approximation. The additional training 
pairs offered no new information to the polynomial approximation. The coefficient v G was 
zero for training pairs using the 3 and 5 factorial designs and a second order polynomial 
approximation. 

2.3 Conclusion 

For a given order of approximation, a good design must use an adequate number of levels 
of the design variables or a poor approximation will be obtained. Likewise, design points 
must be located so that information is available for determining all of the undetermined 
coefficients of the approximating function. In many instances, especially when the region 
of interest is small, a second order polynomial approximation or neural net equivalent will 
be sufficient to build a response surface. A second order approximation requires a design 
containing 3 levels of the design variables. Program DESIGNS has been developed to 
generate a minimum point design which allows all of the coefficients of a second order 
polynomial approximating function to be obtained. This minimum point design can be 
augmented by randomly selected design points or by user selected points. 
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3. Standard Designs 


3.1 Underlying Principle 

When making a polynomial approximation of a function, the number of design levels 
required for each design variable depends upon the order of polynomial approximation 
being used. Consider for example the problem of approximating a function y, a function of 
one design variable. As previously discussed, two levels of the design variable would be 
required to make a linear approximation of the function, three levels of the design variable 
would be required to make a second order approximation, four levels of the design variable 
would be required to make a 3rd order approximation, etc. If y is a function of r design 
variables, a pth order polynomial approximation,^, requires designs at p+1 levels in each 
design variable. 

In response surface methodology, the term factor is used for design variable. A factorial 
design or factorial experiment is a design in which one uses each of the possible 
combinations of the levels of each factor. If m is the number of level of each factor and r 
is the number of factors, then the design would be referred to as a m r factorial experiment . 
Table 3.1 gives the number of designs in various factorial experiments. 
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Table 3.1. Number of designs in a full factorial design 


m = level 
r= factor 

2 

3 

4 i 

2 

4 

9 

16 

3 

8 

27 

64 

4 

16 

81 

256 

10 

1024 

59049 

1.05E06 


One can see that even for a small number of factors, complete factorial experiments become 
impractical if designs are computationally or experimentally expensive to obtain. One then 
is forced to use some sub-set of the factorial design or alternate designs containing requiring 
fewer design points. Concepts from statistics are normally used in selecting a sub-set of the 
factorial design or in developing alternate designs. Thus statistical concepts are reviewed. 

3.2 Statistical Concepts 

When making an approximation, y, of a function, y, most approaches used to select design 
points for a design consider that 

1. polynomial approximations are employed and 

2. the value of the function, y it determined at the designs, {x} s> contains some error, 6;. 
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A measure of the error at point i is the variance of the error, var(e;) = a 2 where 


n 


O 


2 



i*l 




( 30 ) 


where 

H is the true mean of all possible observations of y ; and 
n is the number of observations made. 

In experimental investigations, is experimental error. When making approximations to 
analytical functions, e ; is zero and the variance of the error at point i is zero. Often 
approximations are made to a function whose values must be obtained from some numerical 
algorithm such as the finite element method or finite difference method. Values of y { 
obtained from such algorithms depend on control parameters which dictate the level of 
accuracy of the solution. For example, if y was a stress determined from a finite element 
analysis, then y could depend on a control parameter which specifies the coarseness of the 
finite element idealization. In this case, different values of y s would be obtained for the ith 
design for different values of the control parameters and could be thought of as a 
numerical error. 

It would be an interesting study to select designs such that approximations developed are 
insensitive to numerical errors such as finite element idealization error. However, the 
problem at hand is to find a good approximation to an analytical function or a good 
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approximation for output from a deterministic model. For the problem at hand, for a given 
design, x i5 one obtains the same functional value, yj no matter how many times the function 
is evaluated. Thus, the problems considered in this report contain no numerical error. 
However, as all known algorithms with one exception [21] consider that there is some 
experimental or numerical error, this section now further examines this case. 

Errors in the value of y ( used to build an approximation affect the estimation of the 
undetermined coefficients, bj, in the polynomial approximation and thus affect ft, the values 
of y ; predicted by the approximation. A measure of the error in bj resulting from errors in 
y ; is the variance of bj. For example, consider that y t is obtained from a finite element 
analysis and that a pth order polynomial approximation is employed. The undetermined 
coefficients in that approximations, bj, can be determined from Equation (12). If a number 
of approximations were now made with finite element results, obtained using different 
idealizations, the coefficient bj for these approximations would be different. The variance 
of bj is a measure of how much the b’s change for these different approximations. In like 
form, the different approximations yield different ft and the variance of ft is a measure of 
how much the ft values change from approximation to approximation. 

From a numerical standpoint, it is desirable to have approximations that are not highly 
sensitive to the error Approximations are insensitive to the error, e,, if the variance of 
bj and the variance of ft is small. Most design selection algorithms currently in use attempt 
in some way to keep these variances small. 
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The variance of bj is the j,j term of the variance-covariance matrix cov b where (see 
Equation 3.11 of [3] or Equation 2.8 of [2]) 


[covb^oHmZ ]) 1 


(31) 


and the variance of ft is given by (see Equation 2.11 of [2]) 

va r y i =o 2 [Z i } , ([Z] t mr'{Z i ) (32) 

where {Zy is the lxp vector whose elements correspond to the elements of a row of matrix 

[Z]. 

Notice that these variance involve the matrix [H] where 

[H]=m‘[z]y l (33) 

Design selection affects [Z], which from Equation (33) affects [H], which in turn affects the 
variances of bj and ft. Many design point selection algorithms attempt to select designs 
which give an [H] matrix which will keep the variances of bj and ft small. 

3.3 Orthogonal Designs 

The associated undetermined coefficients of a polynomial approximation function can be 
found from Equation (12). The solution for these coefficients involve the matrix [Z] (see 
Equations (9) and (10)). Let {Z} be the ith column of matrix [Z]. A design is said to be 
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orthogonal if the columns of the [Z] matrix are orthogonal, i.e. {ZjJ'IZj} =0, i^j. There are 
interesting properties of orthogonal designs which have prompted there use. Thus 
orthogonal designs will now to presented in some detail. 

3.3.1 Scaling 

The discussion of orthogonality is simplified by working with scaled variables. Consider that 
the approximation in question involves k unsealed design variables Xj and contains N design 
points. Instead of working with x j5 the variables will be scaled. Let x iu be the uth level of 
unsealed variable i and x iu be the scaled level. The desired scaling is 


N 


E 


x « =0 ’ ,=1 >* 


This scaling can be accomplished by having 





(35) 


(36) 


where 


and 


x=the average of the levels of x { 


(37) 
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With this scaling, N experimental design points of the orthogonal design give 


[Z\'[Z]=m 

N 


where [I] is the identity matrix. 


3.3. 1.1 Example of Scaled Designs: 

Consider a 2 factorial design with levels of 4 and -4. For that design 


and 


ij=0, x 2 =0 


or Si=S2=4 

From Equation (3), the levels of the scaled variables are 


* J ‘*~° 

4 


or the levels of the scaled variables are 1 and -1. 



3.3.2 Bias 


Assume that the polynomial approximating function is inadequate. The coefficients of that 
polynomial can be determined from Equation (12). Let {bj} be the coefficients thus 
obtained and let [Z t ] be the corresponding [Z] matrix. Then from Equation (12) 




( 44 ) 


Assume that the function being approximated can be expressed as 

{Y}=[Z\{b} 


( 45 ) 


where 


{b} = 



[Z\=[ [Z t ] 


[Z 2 ] ] 


( 46 ) 


Entering Equations (40), (45), and (46) into Equation (44) gives 


A 


t-Imizjw,] im 

N 


Pl> 


> 


( 47 ) 


Entering Equation (39) into Equation (47) gives 


( 48 ) 


or 


26 


( 49 ) 


[b l } = {b 1 \ +h Z t y[Z 2 ][b 2 } -1M+M {b 2 } 
N 


where [A] is called the alias matrix . One can see in Equation (49) that the coefficients {bj} 
will only be correct estimates of {bj if the columns of [ZJ are orthogonal to the columns 
of [ZJ. Special situations where this orthogonality occurs are next discussed. 

3.3.2. 1 A bias example— linear approximating polynomial but the exact function contains 
linear terms and cross-product terms: 

Consider a linear approximating polynomial 

(5#) 
i- 1 


where the exact function is 

H + E^ + EEV^ 

i=l i»l /»i 


(51) 


where by are the undetermined coefficients associated with the cross-product terms. For this 
problem, a full 2 k factorial design gives that the columns of [ZJ are orthogonal to the 
columns of [Z 2 ] and thus 
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( 52 ) 


33.2.2 A bias example-linear approximating function but the exact function is a complete 
quadratic polynomial: 


Consider a linear approximating polynomial 


y=k + Y, fa 


(53) 


where the exact function is a complete second order polynomial thus 

y-VE + E Vf + E E 

i=l i=l i«l y»i 


(54) 


Assume again that a full 2 k factorial design is used. For this problem the alias matrix is 
such that one obtains 


WE 


i-1 


bj=b p j=ljc 


(55) 


Thus only 6 0 is biased with the other coefficients unbiased or uncorrelated. 
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3.3.3 Orthogonal Designs for Linear Approximations 

For a problem with r design variables, a full 2 r factorial design is an orthogonal design if the 
approximating function is a first order polynomial. There are several advantages in using 
such an orthogonal design when the approximating function is assumed to be linear. These 
advantages are: 

1. The solution for the coefficients of the polynomial approximation require a matrix 
inverse (see Equation (12)). However, when the design is an orthogonal design, that inverse 
is very easily obtained using Equation (40). Thus there is a small computational advantage 
in using an orthogonal design. 

2. Examples 3.3.2.1 and 3.3.2.2 indicate that under certain conditions, the coefficients 
obtained using an orthogonal design are unbiased. Obtaining unbiased coefficients is 
probably more important in developing response surface from experimental results than 
when developing response surfaces when results are from a deterministic model. With 
experimental studies, it may be important to ascertain the unbiased values of the linear 
coefficients. For the deterministic model however, one is looking for an 
approximating function which gives a good approximation throughout a region of interest. 
Whether the coefficients of the polynomial approximation are biased or unbiased is of little 
concern. 

3. It can be proven that for linear polynomial approximations, an orthogonal design gives 
the minimum variance of the coefficients (see page 109 of [3]). It is important when 
modeling experimental results to obtain a model that is not overly sensitive to experimental 
error and thus there is an advantage in having a minimum variance of the coefficients. 
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However, for response surfaces of a deterministic model, variance of the coefficients is not 
relevant. 

3.3.4 Orthogonal Designs for 2nd Order Polynomial Approximations 
It is not possible to find an orthogonal design when using a second order polynomial 
approximating function of the form of Equation (8) (see page 107 of [2]). However, an 
orthogonal design can be found if one uses as the approximating function a second order 
orthogonal polynomial (page 130 of [3]) 

y= b o + E V + E V*< 2 -*?) + E E (56) 

i=l i«l i«l y*i 


where 



s 



N 


(57) 


and where 

N=the number of design points and 
x. =x, for each of the design points. 


(58) 


The use of an orthogonal design still gives the small computational advantage that the 
inverse shown in Equation (12) is an inverse of a diagonal matrix. However, when using 
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second order approximations, it is not clear under what conditions one obtains unbiased 
coefficients. Also it can not be proven that orthogonal designs any longer give a minimum 
variance of the coefficients. Thus most of the reasons for using orthogonal designs found 
for linear approximations are not present when using second order approximations. 

3.3.5 General Discussion of Orthogonal Designs 

Orthogonal designs offer a small computational advantage that the matrix inverse required 
in solving for the coefficients of the polynomial approximating function is an inverse of a 
diagonal matrix. When approximating a deterministic model, properties of orthogonal 
designs which minimize the variance of the coefficients and which give unbiased coefficients 
are unimportant. For this case, the use of orthogonal designs can only be justified by how 
well they perform on test problems. Such test problems are presented later in this report. 

3.4 Central Composite Designs-Designs for Fitting Second Order Models 

It was shown in Section 2 that at least 3 levels of the design variables are required if one 
is to make a second order approximation. A workable alternative to using a 3 k factorial 
design is a class of designs called the central composite design . These types of designs are 
widely used by workers applying second order response surface techniques [3]. 

3.4.1 Format of the central composite design 

The central composite design is a design composed of the 2 k factorial design augmented by 
additional points. The augmented design points are as follows: 
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*1 *2 *3 • 

0 0 0 . 

-a 0 0 . 

a 0 0 . 

0 -a 0 . 

0 a 0 . 

0 0 0 . 

0 0 0 . 

Figure 7 shows a central composite design for k=3. The value of a and the number of 
design points at the center of the design are varied to meet certain conditions. In the 
following, those conditions are chosen assuming that the approximating polynomial function 
is given by Equation (56). 

3.4.1. 1 Single center point rotatable second order experimental designs: 

A design is said to be rotatable when the variance of the estimated response— that is, the 
variance of y, which in general is a function of position in the design space, is instead only 
a function of the distance from the center of the design and not on the direction. In other 
words, a rotatable design is one for which the quality of the estimator y is the same for two 
points that are the same distance from the center of the design [3]. It is possible to develop 
central composite designs which have a single center point. The value of a which will yield 
these rotatable second order designs are given in Table 3.2. 


0 

0 

0 

0 

0 

-a 

a 


(59) 
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Table 3.2. Value of a for single center point rotatable central composite designs 


k 

a 

2 

1.414 

3 

1.682 

4 

2.000 

5 

2.378 

5 (1/2 rep) 

2.000 

6 

2.828 

6 (1/2 rep) 

2.378 

7 

3.364 

7 (1/2 rep) 

2.828 

8 

4.000 

8 (1/2 rep) 

3.364 


Note in Table 3.2 that a rotatable second order experimental design can be obtained with 
a fractional factorial design augmented with additional design points as well as with a 
augmented full factorial design. 


3.4.1.2 Multiple center point rotatable uniform precision designs: 

In general, the variance of y varies with distance from the center of the design. However, 
by varying the number of center points, N, the variance at a distance of unity from the 
center can be made approximately equal to the variance at the center of the design. Such 
designs are referred to as uniform precision designs . The uniform precision design is based 
on the philosophy that in the central region of the design space there should be uniform 
importance as far as the variance of response is concerned, as opposed to, for example, a 
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situation in which the variance is low in the center of the design but increases drastically as 
one moves away from the design center [3]. The number of center points, m, and the value 
of a can be varied so as to obtain a rotatable uniform precision designs . Table 3.3 gives 
those values. 


Table 3.3. Values of m and a for multiple center point rotatable uniform precision designs 


k 

m 

a 

2 

5 

1.414 

3 

6 

1.682 

4 

7 

2.000 

5 

10 

2.378 

5 (1/2 rep) 

6 

2.000 

6 

15 

2.828 

6 (1/2 rep) 

9 

2.378 

7 (1/2 rep) 

14 

2.828 

8 (1/2 rep) 

20 

3.364 


3.4.1.3 Single center point orthogonal central composite designs: 

An orthogonal central composite design can be developed where [Z] l [Z] is diagonal. To 
obtain a design of this type a single center point can be used and the a value are taken from 
Table 3.4. 
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Table 3.4. Values of a for single center point orthogonal central composite designs 


k 

a 

2 

1.000 

3 

1.216 

4 

1.414 

5 

1.596 

6 

1.761 

7 

1.910 

8 

2.045 


3.4. 1.4 Rotatable orthogonal designs: 

By varying the number of designs at the design center, m, and by selecting appropriate 
values for a, an orthogonal rotatable central composite design can be obtained. Values of 
m and a for such a design are given in Table 3.5. 
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Table 3.5. The value of m and a for multiple center point orthogonal rotatable central 
composite designs 


k 

m 

a 

2 

8 

1.414 

3 

9 

1.682 

4 

12 

2.000 

5 

17 

2.378 

5 (1/2 rep) 

10 

2.000 

6 

24 

2.828 

6 (1/2 rep) 

15 

2.378 

7 (1/2 rep) 

22 

2.828 

8 (1/2 rep) 

33 

3.364 


3.4.2 Discussion of the central composite design 

Orthogonal central composite designs have been shown to give a variance of response 
comparable to that obtained with a full 3 k factorial design. Thus, their use is justified when 
one has experimental error in the response function. Rotatable and uniform precision 
designs attempt to control the response variance. Thus there use is also justified when one 
has experimental error in the response function. However, when building a response surface 
for a deterministic model where there is no experimental error in the response function, 
their use is justified only by how well they perform of trial problems. Likewise, the designs 
were developed for the approximating function of Equation (56). If a different second order 
polynomial approximating function such as in Equation (8) were used or if a neural net was 
used to develop the response surface, then again the justification for the use of the various 
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central composite designs would have to be based on their performance on trial problems. 
Performance of various central composite designs on trial problems is next reported. 


3.4.3 Example - Fox’s Banana Function 
Fox investigated in Reference [16] a function 

y=10xi~20x 2 xf+10x^+xf-2x l +5 

which has banana shaped contours as seen in Figure 8. The region of interest to be 
considered is (-1.5<x 1 <1.5, *.5<x 2 <2.0). 

A second order polynomial approximation is to be made of this function using an orthogonal 
polynomial approximation as in Equation (56). A two variable orthogonal polynomial 
approximation is of the form 

y=b 0 ^b l x 1 ^b^ 2 +b n (x^-xf)+b n (x^-x^)+b l ^ x x 2 < 61 ) 


where 



N 



tt -1 


N 


( 62 ) 


and where 
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N=the number of design points and 


( 63 ) 


x, =x, at the design points 

In the first phase of this example, Fox’s function was approximated using the second order 
orthogonal polynomial of Equation (61). The designs used in making the approximation 
were 

1. a full 5 2 factorial design, 

2. a full 3 2 factorial design, 

3. single center point rotatable central composite design, 

4. multiple center point rotatable uniform precision central composite design, 

5. single center point orthogonal central composite design, 

6. multiple center point rotatable orthogonal central composite design, 

7. minimum point design from program DESIGNS, 

8-10. minimum point design from program DESIGNS augmented by additional randomly 

selected design points, and 

11. nine randomly selected design points. 

Once an approximation was obtained, the approximate function was evaluated at a 31 x 31 
grid of points over the region of interest. The approximate function values at these 961 
points were used to develop the error parameter v G from Equation (6). Because there are 
a differing number of functional evaluations required for each of the sundry designs tested, 
a comparison of the designs based on v G is misleading. For example, the full 5 2 factorial 


38 



design has 25 design points each requiring a functional evaluation where as the multiple 
center point rotatable orthogonal central composite design has but 16 design points 
requiring 9 functional evaluations (in the following it is assumed that the function being 
approximated has no experimental or numerical error and thus the 8 design points at the 
design center require but one functional evaluation). Thus a comparison of performance 
based only on quality of fit is not a fair comparison. The 5 2 factorial might do a better job 
of approximating a function but the computational cost of the 25-9 = 16 extra functional 
evaluations might make it a less desirable design. 



The efficiency of all the designs was compared to design 1, the 5 2 factorial design. Table 
3.6 gives, for each design tested, the number of design points, N; for central composite 
designs, the number of design points at the center of the design, m; the number of 
functional evaluations required, T; the value of v; the value of v G ; and the value of Ej. 
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Table 3.6. Performance of various designs on Fox’s Banana Function, orthogonal 
polynomial approximating function, -1.5<x 1 <1.5, -.5<x 2 <2.0 


Design 

N 

m 

E9 

V 

V G 

mm\ 

5 2 factorial design 

25 

... 

25 

70.76 

78.92 

1.00 

3 2 factorial design 

9 

««• 

9 

64.07 

102.46 

.47 

single center point rotatable 
central composite design 

9 

■ 

9 

54.36 

77.34 

.35 

multiple center point rotatable 
uniform precision central 
composite design 

13 

5 

9 

53.08 

77.34 

.35 

single center point orthogonal 
central composite design 

9 

1 

9 

64.07 

102.46 

.47 

multiple center point rotatable 
orthogonal central composite 
design 

16 

8 

9 

51.62 

77.34 

.35 

minimum point design from 
program DESIGNS 

6 

... 

6 

0 

162.62 

.49 

minimum point design from 
program DESIGNS augmented 
by 2 randomly selected design 
points 

8 

... 

8 

43.27 

105.16 

.43 

minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 

9 

... 

9 

53.53 

88.63 

.40 

minimum point design from 
program DESIGNS augmented 
by 4 randomly selected design 
points 

10 

••• 

10 

53.05 

86.44 

.44 

random-9 points 

9 

m 


21.05 

460.96 

2.10 
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Several items can be noted in Table 3.6: 

1. The design composed of 9 randomly selected design points did poorly. Even though the 
design points were chosen randomly, it turned out that the design points were not well 
scattered in the design space but were heavily concentrated in one quadrant of the design 
space. The polynomial approximation fitted the function well at the design points but poorly 
over the region of interest. 

2. The value of v G was approximately the same for the single center point rotatable central 
composite design, the multiple center point rotatable uniform precision central composite 
design, and the multiple center point rotatable orthogonal central composite design. These 
three designs differ only in the number of design points at the center of the design space. 
These designs have 1, 5, and 8 designs at the center, respectively. The effect of putting 
more designs at the center is to translate the response surface toward the center response. 
For this problem, however, the actual and approximated response were very close at the 
design center point, even for only 1 design point at the center. Thus, adding more design 
points at the design center did little to translate the response surface and thus did not 
material effect the value of v G . 

3. The eleven designs of Table 3.5 were next used to build an approximation using the 
standard second order polynomial approximation of Equation (8) instead of the orthogonal 
polynomial approximation of Equation (61). Results identical to those of Table 3.5 were 
found. The type of approximating polynomial may effect variances but does not affect 
quality of fit at the design points or over the region of interest. For those problems were 
there is no experimental or numerical error associated with functional evaluations, one is 
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not interested in variance. Thus, there is little advantage in using the orthogonal polynomial 
approximating functions over a standard second order polynomial function. 

4. Based on efficiency, the single center point rotatable central composite design, the 
rotatable uniform precision central composite design, and the rotatable orthogonal central 
composite design performed the best but none of the designs gave a good approximation 
over the region of interest. Over a small region of interest, one could expect that a second 
order polynomial approximation could well approximate the given function. Obviously, here 
the region of interest is too large for a second order approximation to be a good one. Thus 
a smaller region of interest was chosen, -.5 <x 1 ,.5, -.5 <x 2 < .5. Table 3.7 compares the eleven 
designs using this region of interest. Notice that over this smaller region of interest, all the 
designs gave a much better approximation to the function. 

5. For the smaller region of interest, based on efficiency, the 3 2 factorial design, the single 
center point orthogonal central composite design, and the augmented minimum point 
designs performed the best. Obviously, the optimum choice of design is problem dependent. 
However, all designs except the randomly selected design performed much better than the 
5 2 factorial design. 
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Table 3.7. Performance of various designs on Fox’s Banana Function, orthogonal 
polynomial approximating function, -.5<Xj<.5, -.5<x 2 <.5 


Design 

N 

m 

T 

V 

v G 

E i 

5 2 factorial design 

25 

««« 

25 

11.16 

8.57 

1.00 

3 2 factorial design 

9 

... 

9 

13.27 

10.95 

.46 

single center point rotatable 
central composite design 

9 

i 

9 

6.58 

14.74 

.62 

multiple center point rotatable 
uniform precision central 
composite design 

13 

5 

9 

5.88 

14.74 

.62 

single center point orthogonal 
central composite design 

9 

1 

9 

13.27 

10.95 

.46 

multiple center point rotatable 
orthogonal central composite 
design 

16 

8 

9 

5.47 

14.74 

.62 

minimum point design from 
program DESIGNS 

6 

■ 

6 

0 

18.66 

.52 

minimum point design from 
program DESIGNS augmented 
by 2 randomly selected design 
points 

8 

I 

8 

5.74 

11.82 

.44 

minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 

9 

| 

9 

6.45 

10.53 

.44 

minimum point design from 
program DESIGNS augmented 
by 4 randomly selected design 
points 

10 

| 

10 

6.33 

10.29 

.48 

random-9 points 

9 

m 

9 

2.42 

47.22 

1.98 
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3.4.4 Conclusion 

Second order polynomial approximations or neural net equivalents are often adequate for 
building response surfaces, especially if the region of interest is small. Central composite 
designs are convenient for building the second order approximations. They provide the 
necessary information for determining all of the coefficients of the approximating polynomial 
and give a good distribution of points in the design space. The approximating function can 
be made to closely fit the exact function at the design center by using multiple center points. 
When modeling deterministic systems, each functional evaluation at the design center yields 
the same function value. Thus, for deterministic models, only one functional evaluation 
need be performed at the center point even when multiple center points are used. Table 
3.8 gives information relevant to central composite designs for various number of design 
variables, k. Central composite designs give over-determined second order polynomial 
approximations. In other words, there are more design points in the design than there are 
undetermined coefficients in a second order polynomial approximation. Table 3.8 also gives 
the percentage that the approximation is over-determined. Previous studies [19] have 
indicated that designs which give approximations that are around 20-50% over-determined 
tend to be efficient designs. One can see that the central composite designs are reasonable 
for k<6. For larger k values, too many design points are being used by the central 
composite designs. For k>5, an augmented minimum point design is a better choice. 
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Table 3.8. Information relevant to central composite designs for various number of design 
variables 


Number of design 
variables, k 

Number of 
coefficients in a 
2nd order 
polynomial 
approximation 

Number of 
functional 
evaluations 
required with a 
central composite 
design 

% over-determined 

1 

3 

4 

33 

2 

6 

8 

33 

3 

10 

14 

40 

4 

15 

24 

60 

5 

21 

42 

50 

6 

28 

76 

171 

7 

36 

142 

294 

8 

45 

272 

504 
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4. Optimality Criteria 


4.1 D. A. E. G. and V Optimality Criteria 

It was pointed out in Section 3 that even for a small number of factors, a complete factorial 
experiment become impractical if functional evaluations are computationally or 
experimentally expensive to obtain and thus one is forced to use some sub-set of the 
factorial design or an alternate design requiring fewer experiments. Section 3 shows that 
the variances of the coefficients of a polynomial approximation and the variance of the 
predicted response involve the matrix [H] given in Equation (33) and repeated here: 

[HMimz\r' (65) 

Schoofs [22] lists five criteria for selecting a sub-set of the factorial designs. These criteria 
involve the matrix [H]. The criteria, referred to as optimality criteria, attempt to make [H] 
minimal. However, "the minimum of a matrix is not a well defined concept and a number 
of operational criteria have been developed" [22]. The optimality criteria for selecting a 
subset of a full factorial design can be based on selecting the subset satisfying the following 
criteria: 

1. D-optimality, which is achieved if the determinant of [H] is minimal which in term gives 
that the product of the eigenvalues of [H] is minimal. 

2. A-optimality, which is achieved if the trace of [H] is minimal which in term gives that 
the sum of the eigenvalues of [H] is minimal. 

3. E-optimality, which is achieved if the largest eigenvalue of [H] is minimal. 


46 




4. G-optimality, which is achieved if the maximum over all candidate points of the 
estimated response variance is minimal. 

5. V-optimality, which is achieved if the estimated response variance, averaged over all 
candidate points is minimal. 

4.1.1 Criteria Applied to a One Dimensional Example 

An example is considered here to compare the performance of the 5 optimality criteria. 
The following test function of one variable was considered: 

y=2+x+sin[— (x+1)], -l^c^l (66) 

2 


This function was approximated with polynomials of order 1-4. The approximations shown 
in Figure 9 were developed using 13 designs, uniformly spaced in the region of interest. 
These approximations were then used to generate the functional values at 61 uniformly 
spaced points in the region of interest which were used to plot the curves of Figure 9. 

Further approximations of Equation (66) were developed using various number of design 
points, n. The designs selected were 

1. uniformly spaced design points, n= 5,7,9, 11, 13; 

2. randomly selected design points, n= 5,7,8,11,13; 

3. an n member subset of the 13 uniformly spaced design points, n= 5,7,9, 11. 
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Under item 3, the subset of design points was chosen using: 

1. D-optimality, 

2. A-optimality, 

3. E-optimality, 

4. G-optimality, and 

5. V-optimality. 

A FORTRAN program was written to perform the investigation under item 3. The 
demanding part of the programming was to identify all the possible subsets from the set of 
thirteen design points. After developing a procedure to identify all combinations, each 
subset was used to build the [H] matrix. The "optimal" [H] matrix was then determined 
using the five optimality criteria. The coefficient v G was then computed for the optimal 
subset. Figures 10-13 show the value of v G for the D, A, E, and G optimality criteria when 
a first, second, third, and fourth order approximation is being made, respectively, versus the 
number of design points specified in the subset . Also shown in those figures is the value 
of v G for designs consisting of design points uniformly spaced in the region of interest. 

It was found that for all subsets of size r from a design point set of size n that the estimated 
response variance, averaged over all candidate points, was invariant. This finding 
undoubtedly could also be proven theoretically but such a proof was not attempted. From 
this example, one can conclude that the V optima lity criteria, which employees the estimated 
average response variance, is not a viable criteria for sele ct i ng a subset of design points from 


48 



a given set . From Figures 10-13, one can see that in most cases there is little difference in 
the performance of the various optimality criteria with criteria D and G performing slightly 
better than the other two criteria. As can be seen in Figure 12, on one occasion (when 
using a third order polynomial approximation and when selecting a subset of 5 design points 
from the 13 design point set) the G optimality criteria performed poorly while the D criteria 
did not. Thus, this example indicates that the D optimality criteria mav be the criteria_of 
choice . There is a further advantage in using the D optimality criteria. The requirement 
that the determinant of [H] is minimal is equivalent to a requirement that the determinant 
of [G] is maximal where 

[GMZftZ] (67) 

Thus the D optimality criteria insures that the procedure for determining polynomial 
coefficients in Equation (12) will be well defined. In other words. Equation (12) uses the 
inverse of [G]. The D optimality criteria guarantees that [G] is not singular. 

One can see in Figures 10-13 that, in most cases, all the optimality criteria performed worst 
than the uniformly spaced design case. This example indicates that a design picked Vising 
an optimality criteria may be no better than a de s ign of the same size in which the design 
points are uniformly located in the design space . 

4.2 S and O Optimality Criteria 

The previous optimality criteria involved only the matrix [H] and did not consider the 
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function to be approximated. Thus for a given number of design variables and level of 
approximation, the same designs would be selected no matter what the nature of the 
function to be approximated. Initially it was thought that a superior optimality criteria 
would have to consider the nature of the function. Thus two additional optimality criteria 
were examined: 

1. S-optimality, which is achieved if the average error of approximation at the design points 
is minimal and 

2. Q-optimality, which is achieved if the maximum error of approximation at the design 
points is minimal. 

Here 

( 68 ) 

average error of approximation 


and 


maximum error of approximation=tnax (y-y) 2 , i=l,..,r 


(69) 


where r is the size of the subset of design points to be selected. One can see that with the 
S and Q optimality criteria, the function to be approximated effects the design points 
selected. 
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4.2.2 Criteria Applied to a One Dimensional Example 

The one dimensional example problem of Section 4.1.1 was then re-examined. Figures 14- 
17 show values of v G using the S and Q optimality criteria and using a first, second, third, 
and fourth order polynomial approximation, respectively, versus size of the subset of design 
points. Also shown in these figures are results for uniformly spaced design points. One can 
see in these figures that terrible approximations were obtained with these criteria when only 
small subsets of design points were selected from the original set. Figures 18-20 indicate 
why such bad approximations are obtained with these two criteria. 

Figure 18 depicts results obtained by having eleven designs points selected, using the Q 
optimality criteria, from a set of 13 design points. The Q optimality criteria finds an 
approximation such that the maximum error of the approximation over eleven design points 
is minimal. One can see in Figure 18 that the approximating function did indeed well fit 
the exact function at the 11 design points selected. However, the approximating function 
did a poor job of approximation at the ends of the region of interest and thus would not 
yield a low value of v G . Figure 19 is similar to Figure 18 except that this figure depicts 
results obtained by having 7 design points selected from the set of 13 design points. One 
can see that for the optimum design selected, there is an almost perfect approximation at 
the design points selected but over a much larger region the approximation is poor and thus 
a large value of v G would be obtained. In Figure 20, only 5 design points are selected. 
Again at those design points, an almost perfect approximation is obtained but a terrible 
approximation is obtained over a large part of the region of interest and thus a large v G 
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would be obtained. Thus we can conclude that the S and O optimality criteria are not 
operative . 


4.3 An Alternate Approach-Random Selection of Designs 

The effect of randomly picking design points was next considered. Here designs are picked 
in the region of interest using a random number generator. 

4.3.1 Random Selection of Designs Applied to a One Dim ensional Example 
For the one dimensional problem under consideration, first, second, third, and fourth order 
approximations were considered. Design point sets containing 5,7,9,11, and 13 design points 
were developed by randomly picking design points in the region of interest using a random 
number generator. Approximations were developed using the design sets. Results using 
these approximations are compared in Figures 21-24 to results using uniformly spaced design 
points. One can see in these figures that most of the time results from randomly picked 
design points are either as good as or not much worst than results from uniformly spaced 
design points. However, on two occasions, when the number of design points in the design 
set was small, a relatively poor approximation was obtained. Obviously where one is picking 
only a small number of points using a random number generator, there is a chance that a 
bad set of points can be generated and indeed on these two occasion a poor selection of 
points was made. In general however, when more design points are randomly selected, 
those points should be scattered throughout the design space and good approximations 
should be obtained. In conclusion, randomly selecting de sign points mav be a viable method 
of design selection . 
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4.4 Larger Problems 

Consider a problem in two variables and consider that the potential design points will be 
taken from a 6 x 6 grid of points. Let 

r = total number of design points in the set of potential design points, 
c = number of design points in the selected subset of design points, 
nc = the number of different combinations of designs in the subset. 


For the problem at hand, r=36. Subset sizes of c = 15, 20, 25, and 30 are to be considered. 
The number of possible combinations of design points in the subset, nc, is given by 


nc= 


r 1 

(r-c)! cl 


(70) 


Table 4.1 summarizes the number of combinations for this study. 


Table 4.1 Number of combinations of designs in a two variable study 


r 

Total number of design 
points 

c 

Number of point in subset 

nc 

Number of combinations 

36 

15 

5,567,902,560 

36 

20 

7,307,872,110 

36 

25 

600,805,296 

36 

30 

1,947,792 
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One can see that for even small problems, it is infeasible to examine all possible 
combinations of subsets of size N from a given set of design points. Welch [23], instead of 
evaluating all possible N-point designs, developed a "branch and bound” algorithm which 
guarantees global D-optimal designs but which does not generate and evaluate all possible 
designs. However, even here the computing costs are high. Fedorov [24] developed another 
technique which neglects the integer character of the components of the design set and 
obtains a discrete design which is rounded off to an exact design. Reference [22] reports 
that these designs are considered only approximate. The most popular algorithm seems to 
be DETMAX by Mitchell [25]. Quoting reference [22], "The algorithm starts with an initial 
m-point ED (experimental design); the final goal is an optimal N-point ED. During each 
iteration step that candidate point, which results in the largest increase of det(M), is added 
to the design, and subsequently that point, which results in the smallest decrease of det(M), 
is removed from the design. The number m of points in the initial design may be larger or 
smaller than N. If necessary the algorithm first adds (if m<N) or rejects (if m>N) points 
until the number of points in the ED is equal to N. In order to avoid local optima the 
algorithm is able to perform ’excursions’, in which several points are added at one go and 
subsequently the number of points is reduced to N. If the resulting N-point ED has not 
been improved, another excursion will be made from the same initial design. If the 
excursion is successful the resulting ED will be used as starting ED in a further attempt to 
maximize det(M). The algorithm terminates when, after several excursions, no better ED 
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is found. The algorithm generates high quality EDs against relatively low computing costs." 
An attempt is being made to obtain the algorithm DETMAX. 

4.5 Optimality Criteria Based on Minimizing Uncertainty 

Reference [21] considers problems where there is no experimental error. That reference 
uses an optimality criteria based on selecting a design which minimizes the uncertainty in 
the approximating function. That reference was given mixed reviews by a number of leading 
authorities in the field [21] (reviews follow the paper). The formulation is quite theoretical 
and difficult to follow. The formulation seems to have promise but requires additional 
theoretical development before it becomes operative. 

4.6 Conclusion 

There is little rational for using any of the investigated optimality criteria when building 
approximations of functions which contain no experimental error. However, the D- 
optimality criteria can conveniently be used as a heuristic in selecting design points. 

Previous investigations have indicated that approximations should be over-determined. That 
is to say that more training pairs should be used to build an approximations than the 
number of associated undetermined parameters. It has been suggested that a 20-50% over- 
determined system might be reasonable. The program DESIGNS, described in Section 2, 
develops enough designs to exactly determine a quadratic approximation of a given function. 
The D-optimality criteria can be used as a heuristic for selecting design points to 


55 



supplement those generated by DESIGNS. The use of the D-optimality criteria to select 
the supplementary points would guarantee than no singular matrices would be encountered 
in determining the undetermined parameters associated with the polynomial approximation. 
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5. Significance Testing of Coefficients 


5.1 Introduction 

When the training pairs used to build a polynomial response surface contain experimental 
or numerical error, certain coefficients in the polynomial approximation may not be 
significant. In other words, even though one calculates a value for some coefficient, b j5 the 
experimental or numerical error may have such an effect on that coefficient that it could just 
as well be taken as zero as the value calculated. In situations like this, it may be 
advantageous to drop that term from the polynomial approximation and redevelop the 
response surface. Such a procedure is discussed in pages 34-38 of [3] and an automated 
procedure for performing such an operation was developed in [26]. Testing of significance 
involves the t-test which is next described. 

5.2 t-test 

Coefficients of the polynomial approximation are found from Equation (12). The 
determination of those coefficients involve the matrix [H] where 

[H]<[ mz })- 1 (71) 


A number of terms must now be defined: 


mean square error=MSE- 


E 


i-l 


N-m 


(72) 
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standard error coefficient =se i =^MSE H u 


(73) 



(74) 


where 

N = the number of design points and 

m = the number of coefficients in the polynomial approximation. 

In making the test of significance, t ( from Equation (74) is compared to tabulated values of 
t a . The value of t a is taken from a table of "Percentage Points of the Student’s t 
Distribution" [3]. The value taken depends on the level of significance desired. In lieu of 
using tabulated values, t a is often taken as four [26]. If t 4 is less than t a (q < t a ), then that 
coefficient’s importance in approximating the response is deemed to be insignificant and 
therefore may be eliminated from the response function. 

The primary focus of this study was to examine methods of developing good response 
surfaces for deterministic models, i.e. for systems that contain no experimental or numerical 
error. Statistical testing of coefficients presupposes experimental or numerical error and 
thus is not relevant when approximating response which contains no error. However, the 
method was thought to perhaps offer a heuristic for improving the quality of a response 
surface even if experimental or numerical errors are not present. Thus, two examples were 
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examined. Results are next reported. 


5.3 Example 1 - Fox’s Banana Function 

Example 1 again examines Fox’s Banana Function [16]. A complete second order 
polynomial approximation (m = 6) and a complete third order polynomial approximation 
(m = 10) were developed. These approximations were developed using a complete 6 2 
factorial design (N = 36). A t-value, was calculated for each parameter, b„ and compared 
to t a = 4. Parameter that lack significance (t; < t a ) were eliminated. A new approximation 
was then developed using only the significant parameters. The values of v and v G from 
Equations (5) and (6), respectively , were developed for the complete polynomial and for 
the polynomial containing only terms deemed significant. Results are shown in Figures 25 
and 26. On can see in these figures that eliminating coefficients deemed insignificant had 
an adverse effect on the quality of the approximation over the region of interest. 

5.4 Example 2 

The effect of eliminating coefficients deemed insignificant was tested on the function 

T=(4+x 1 ) 3 +sin[— (Xj + 1)] +2+*2 +sin(^-)+7x 2 x 1 (75) 


Again, a complete second order polynomial approximation (m = 6) and a complete third 
order polynomial approximation (m=T0) were developed. These approximations were 
developed using a complete 6 2 factorial design (N=36). A t-value, t j? was calculated for 
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each parameter, bj, and compared to t a — 4. Parameter that lack significance (tj < t a ) were 
eliminated. A new approximation was then developed using only the significant parameters. 
The values of v and v G from Equations (5) and (6), respectively , were developed for the 
complete polynomial and for the polynomial containing only terms deemed significant. 
Results are shown in Figures 27 and 28. On can see in these figures that eliminating 
coefficients deemed insignificant offered no improvement in the quality of the response 
surface. 

5.5 Conclusion 

The applicability of significance testing of polynomial coefficients when modeling 
deterministic systems was considered. Two examples were examined to see if eliminating 
terms of polynomial approximations which were deemed to be insignificant by the t-test 
would improve the quality of the response surfaces developed. Based on these two 
examples, it was concluded that no improvement in the predictive capability of response 
surfaces over regions of interest would be obtained with such a procedure. The relevance 
of significance testing is when modeling systems containing numerical or experimental error. 
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6. Applicability of the Response Surface Technique 


6.1 Introduction 

The following study was performed to ascertain under what circumstances could the 
response surface technique be used to advantage in engineering optimization application. 
In this regard, assume that a quadratic polynomial approximations is to be made of functions 
of n variables. The number of undetermined coefficients in that approximation is: 

number of coefficients^ 


Previous studies [19] have shown that the best approximations are obtained when the 
approximations are over-determined. Thus, the number of functional evaluations required 
to make the approximation is: 

number of functional evaluations = — ( 77 ) 


where S determines the degree that the approximation is over-determined. 

The functional evaluations required to build the approximation are initially performed 
before the start of the optimization process. By using parallel processing, these functional 
evaluations may be less computationally expensive than evaluations made sequentially in a 
direct optimization procedure. The number of required evaluations of Equation (77) is then 
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equivalent to a reduced number of sequential evaluations thus: 


equalivent number functional evaluations^ 


$P(/i+l)(w+2) 

2 


(78) 


where /3 is a coefficient of efficiency associated with parallel processing. 


An optimum solution can be attempted using the response surfaces developed instead of the 
original functions. However, because of the inexact nature of the approximations, a new set 
of response surfaces may have to be developed at the most recent approximate solution and 
another optimal solution attempted. This procedure may have to be repeated a times to 
reach the optimum solution for the original problem. The total number of equivalent 
functional evaluations performed in reaching this optimum is: 


total equivalent functional evaluations = 


gP&(f»+l)(/t+2) 

2 


(79) 


If the solutions was attempted by direct optimization techniques instead of using response 
surfaces, Barthelemy [27] states that a solution can be obtained in most cases using no more 
than i|r first derivative evaluations. If the first derivatives are obtained by finite difference 
formulae, an estimate of the number of functional evaluations required by a direct solution 
procedure is: 

functional evaluations direct methods*ty(n+ 1) 
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If the response surface technique is to be competitive with the direct solution technique, 
then from Equations (4) and (5) one must have: 


2 


( 81 ) 


where y is a convenience factor associated with using response surfaces. In other words, an 
investigator may tolerate more functional evaluations with the response surface technique 
than with the direct solution procedure just for the convenience of using response surfaces. 
Rearranging Equation (81) gives 

[,»!][ (82) 
2y 


Since (n+1) is positive, one obtains 

aP6(/i + 2) Q (83) 

2y 


or 


7t£ 


2»Y 

ap6 


-2 


( 84 ) 
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In review 


oc — number sequential optimizations 
P= parallel processing coefficient 
6 = overdetermined coefficient 

Y = convenience coefficient 
i|r= direct solution coefficient 


Reasonable ranges of the parameters involved are 


a = 1.00-4.00 
p=0.10-1.00 
6 = 1.25-1.75 
Y = 1.00-3.00 

i}» =6.00- 10.0 


( 86 ) 


For an approximate upper bound on the number of design variable that could be 
economical used with the response surface technique take: 


<x = 1.00 

p=o.io 

6 = 1.25 
y=3.00 
f=10.0 


(87) 


giving 


ns498 


( 88 ) 


Under the most unfavorable set of circumstances, that is: 
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( 89 ) 


a =4.00 
p = 1.00 
6 = 1.75 
Y = 1.00 
f=6.00 


one obtains 


n- 0 


(90) 


Thus depending upon the problem, one could use the response surface technique for n-0 
to n=500 variables. Consider the following reasonable set of parameters 


a=3.00 
P=0.50 
6 = 1.25 
y = 1.50 

i|f=8.00 


(91) 


giving 


nsl3 


(92) 


Thus, it is reasonable to assume that the response surface technique could be used for up 
to 10-15 design variables. 
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6.2 Conclusion 

Under the most advantageous circumstances, the response surface technique applied to 
engineering optimization application could be used for up to 500 design variables. Under 
the worst set of circumstances, it is entirely inappropriate. Under normally expected 
circumstances, this technique might be used to advantage for 10-15 design variables. 
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7. Additional Examples 


7.1 Introduction 

The next several examples examine the effect of design selection on the quality of 
approximations. In each case, a second order polynomial approximation is made of a trial 
function. Different number of design variables are considered in each example. Thus, for 
each example different designs are appropriate. In the first example, there are 4 design 
variables. When there are fewer than 6 design variables, central composite designs are a 
possible appropriate choice. Other choices are the 3 k factorial design, the minimum point 
design, the augmented minimum point design, or randomly selected design. All of these 
designs are considered in that example. In the second and third examples, there are 15 and 
20 design variables, respectively. Here, the 3 k factorial design and central composite designs 
contain too many design points to be practical. For these examples, the minimum point 
design, the augmented minimum point design, and the randomly selected design are 
appropriate and are considered. 

7.2 The 35 Bar Truss with 4 Design Variables 

In many response surface applications, the function to be approximated is a relatively 
smooth function of the design variables which can be approximated with a lower order 
polynomial or an artificial neural net with only a few nodes on the hidden layer. A problem 
of this type is shown in Figure 29. In this example, all loads shown in Figure 29 are in kips, 
all members of the lower chord of the truss are assumed to have area, A lf and all members 
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of the upper chord to have area, A 2 , all vertical and diagonal members to have area, A 3 . 
The depth of the truss is H. A response surface is to be constructed for the stress in 
member BC in terms of the design variables, x ; thus 


x r \ IA t , i= 1,3 
x 4 =.09375tf-.4375 


The region of interest is 


2 in 2 <A t < 8 in 2 
6 ft zHz 10 ft 


or in terms of the design variables 

.125^c,.^.5 


(93) 


(94) 


(95) 


A number of designs were used to develop a second order polynomial approximation for the 
stress in member BC. Each approximation was then used to predict stress on a 5 x 5 x 5 
x 5 grid of points. The predicted stress and the actual stress on these NG = 625 grid of 
points were then used to develop v G from Equation (6). The parameter v G is a measure of 
the quality of the approximation over the region of interest. 

The different designs examined required different numbers of functional evaluation. So as 
to get a measure of the quality of fit of the approximation over the region of interest which 
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Table 7.1 The 35 bar truss with 5 design variables, 2nd order polynomial approximation 


— 

Description 

m 

a 

n 

F 

v (%) 

v G (%) 

E . 


3 4 factorial design 

... 

... 

81 

81 

3.34 

2.41 

1.00 


single center point rotatable 
central composite design 

i 

2.000 

25 

25 

0.66 

2.67 

0.34 


multiple center point rotatable 
uniform precision central 
composite design 

■ 

2.000 

31 

25 

0.59 

2.67 

0.34 


single center point orthogonal 
central composite design 

i 

1.414 

25 

25 

1.47 

2.37 

0.30 


multiple center point rotatable 
orthogonal central composite 
design 

12 

2.000 

36 

25 

0.55 

2.67 

0.34 


minimum point design from 
program DESIGNS 

... 

... 

15 

15 

0.00 

3.99 

0.31 

- 

minimum point design from 
program DESIGNS augmented 
by 3 randomly selected design 
points 

... 

... 

18 

18 

0.40 

3.86 

0.36 


minimum point design from 
program DESIGNS augmented 
by 6 randomly selected design 
points 

... 

... 

21 

21 

0.38 

3.91 

0.42 

— 

minimum point design from 
program DESIGNS augmented 
by 9 selected design points 

... 

... 

24 

24 

0.41 

3.77 

0.46 

— 

randomly selected design 

... 

... 

25 

25 

0.00 

824.2 

105 


m= number of design points at the center of the design space 
T = the total number of design points 
F = the number of functional evaluations required 
a = parameter which defines location of certain design points 
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takes into account the number of functional evaluations performed, the efficiency, Ej, from 
Equation (64) was developed for each design. Table 7.1 reports for each design considered, 
the efficiency, Ej, as well as other relevant information. 

One can see in Table 7.1 that all the designs considered, except the randomly selected 
design, gave a good approximation over the region of interest. Randomly selected designs, 
which often work well, can sometimes suffer from the problem that the coefficient matrix 
used to solve for the approximation’s associated parameters is poorly conditioned or that 
the design points selected are not well scattered throughout the design space. In either case, 
they can yield a poor approximation over the region of interest as in this example. 

The 3 4 factorial design well approximated the trial function. However, because it uses so 
many design points its efficiency measure is poor and thus is not a design of choice. The 
single center point orthogonal central composite design and the minimum point design from 
program DESIGNS performed the best, based of their efficiency. However, excluding the 
randomly selected design and the 3 4 factorial design, all of the designs considered gave a low 
value of v G and had approximately the same value of efficiency. 

Under normal circumstances, information is not available to calculate v G and one must use 
the parameter v as a measure of the quality of fit over the region of interest. However, the 
parameter v is only a measure of quality of fit over the region of interest if the 
approximation is over-determined. Thus, under normal circumstances one would not want 
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to use the minimum point design. This example indicates, that for problems of the size of 
this example, that any of the central composite designs or the augmented minimum point 
designs would be appropriate. 

7.2 The 35 bar truss with 15 design variables 

This example again considers the 35 bar truss of Figure 29. In this example, H is 10 ft., the 
areas of the 14 bars of the top and bottom chords are Aj, i = 1,14, and the area of the 
vertical and diagonal members is A 15 . The design variables of the problem are taken as 

x r l IA P i=l,15 ( 96 > 


The region of interest is 


2 in 1 <: A i <; 8 in 2 


(97) 


or in terms of the design variables 

. 125^.5 ( 98 > 

Response surfaces were developed for the stress in member BC using a 2nd order 
polynomial approximation. The approximation were developed using various designs. To 
test the quality of the approximations over the region of interest, the function and the 
approximations were evaluated at NG=500 randomly selected test points over the region 
of interest. That information was then used to calculate v G from Equation (6). The random 
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number generator used to develop design points uses, in generating its numbers, an initial 
seed parameter, IFLAG. A different value of IFLAG was used to generate the 500 test 
points than was used to generate random points in the randomly selected designs or in the 
augmented minimum point designs. Thus, the test set of points does not duplicate any of 
the design points in the designs considered. Results of this investigation are reported in 

Table 7.2. 

One will notice in Table 7.2 that only minimum point designs, augmented minimum point 
designs, and randomly selected designs are considered. A 3 15 factorial design contains over 
14 million design points. Thus, the use of the 3 15 factorial design is out of the question. For 
a problem in k design variables, the central composite design uses a 2 k factorial design 
augmented by 2k +1 additional design points. Thus, such a single center point central 
composite design for this problem contains 32,799 design points. Here again, such a design 
is impractical. One can develop a central composite design by augmenting only a fraction 
of the 2 k factorial design. For this problem, a single center point central composite design 
using only a 1/4 fraction of the 2 15 factorial design would contain 8,223 design points which 
is still an impractical design. Thus, for problems of the size of this example, only the 
minimum point designs, augmented minimum point designs, and randomly selected designs 
are of reasonable size. 

We can see in Table 7.2 that all of the designs with the exception of the "randomly selected- 
-exactly determined design" did a good job of approximating truss behavior. A singular 
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matrix was encountered in Equation (10) for the randomly selected— exactly determined 
design. With completely randomly selected designs, there is always the possibility of having 
a poorly conditioned coefficient matrix [Z] in Equation (10) and indeed this occurred in this 
problem. However, there was no problem with matrix conditioning using randomly selected 
over-determined designs. 


Table 7.2 The 35 bar truss with 15 design variables, 2nd order polynomial approximation 


Description 

F 

v % 

v G % 

E . 

minimum point design from 
program DESIGN- 
exactly determined 

136 

0 

1.263 

1.0 

augmented minimum point 
design— 20% over- 
determined 

163 

0.083 

0.294 

0.28 

augmented minimum point 
design-40% over- 
determined 

190 

0.087 

0.060 

0.07 

random selection-exactly 
determined 

136 

* 

* ^ 

* 

random selection— 20% 
over-determined 

163 

0.003 

0.029 

0.03 

random selection— 40% 
over-determined 

190 

0.003 

0.010 

0.01 


* singular coefficient matrix 


The efficiency parameter, Ej, is calculated in Table 7.2 but it is rather a meaningless 
parameter for this problem because all the designs so well fit the exact function. In real life 
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situations, one usually does not have available information for calculating v G . Thus, the 
parameter v or like term must be used as a measure of the quality of the approximation. 
The parameter v is not a meaningful measure of the quality of fit over a region of interest 
unless the system is over-determined. Thus for this example, the design of choice would be 
either the 20% over-determined minimum point design or the 20% over-determined 
randomly selected design. 

7.3 Analytical function--20 design variables 

This example considers a problem with even more design variables. The function tested is: 

20 20 20 20 20 

(99) 

f=l i* 1 j*i 1*1 j=i 


A second order polynomial function was used to build the response surface approximating 
this function. The polynomial approximating function had 231 undetermined coefficients. 
Because of the large size of this problem, factorial designs and central composite designs 
are not appropriate. A minimum point design, augmented minimum point designs, and 
randomly selected designs were considered. Values of the test function and approximate 
function were evaluated at NG = 1000 randomly selected points and the parameter v G was 
developed using this information. The measure of efficiency of the designs examined along 
with other relevant information is given in Table 7.3. 
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Table 7.3 Analytical function with 20 design variables, 2nd order polynomial approximation 


Description 

F 

v % 

v G % 

E i 

minimum point design from 
program DESIGN- 
exactly determined 

231 

0 

88.93 

1.0 

augmented minimum point 
design-20% over- 
determined 

277 

5.83 

49.82 

0.67 

augmented minimum point 
design--40% over- 
determined 

323 

9.58 

18.03 

0.28 

random selection-exactly 
determined 

231 

* 

* 

* 

random selection-20% 
over-determined 

277 

0.61 

7.21 

0.10 

random selection-40% 
over-determined 

323 

0.46 

1.20 

0.02 


’ poorly conditioned coefficient matrix 


Just as in Example 7.2, a exactly determined randomly selected design gave a poorly 
conditioned coefficient matrix. These examples indicate that randomly selected exactly 
determined designs should be avoided. The 40% over-determined randomly selected design 
did an excellent job of modeling the test function and was the most efficient design 
considered. It seems that on problems with a large number of design variables that 
randomly selected over-determined designs should be expected to work well. 
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7.4 Conclusion 

The examples of this section have shown that design selection depends on the number of 
design variables. If the number of design variables is less than 6, appropriate designs are: 

1. augmented minimum point designs 

2. central composite designs 

3. over-determined randomly selected designs. 

When there are more than 6 design variables, the central composite designs contain too 
many design point for consideration. For more than 6 design variables, appropriate designs 
are then 

1. augmented minimum point designs 

2. over-determined randomly selected designs. 

The example examined indicate that in all cases, over-determined designs should be used. 
They the most efficient designs. Also, when a design is over-determined the coefficient v 
can be used as a measure of the quality of the approximation over a region of interest. 
Being able to use v as a measure of the quality of fit over the region of interest is very 
im portant because, in general, information is not available to determined the parameter v G . 
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8. Augmented Minimum Point Designs 


8.1 Introduction 

Design selection in the literature concentrates of linear or quadratic response surfaces. This 
study has also concentrated on quadratic approximations for several reasons: 

1. linear approximations, in most instances, will be inadequate to model functions of 
interest, 

2. for many problems, a 2nd order approximation will be adequate to model response 
especially if the region of interest is limited, 

3. there is a scarcity of literature which address design selection for cubic or higher order 
polynomial approximations, and 

4. in optimization process using response surfaces, for moderate size problems, it is more 
computationally efficient to perform a sequence of quadratic approximations than one cubic 
or higher order approximation. This fact is next discussed. 


The number of terms in a second order polynomial in n design variables is 


number terms quadratic =(n+ 


n, w(n+l) 
2 


( 100 ) 


The number of terms in a 3rd order polynomial in n design variables is 


77 



n\ 

6(n-3)! 


(101) 


number terms cubic=l+—n(n+l)+ 

2 


Table 8.1 gives, for various number of design variables, the number of terms in a 2nd order 
and 3rd order polynomial and their ratio. 


Table 8.1 Number of terms in a 2nd and 3rd order polynomial and their ratio 


number of design 
variables, n 

number of terms 
in quadratic 

number of terms 
in cubic 

cubic/quadratic 

3 

10 

20 

2 

6 

28 

84 

3 

9 

55 

220 

4 

12 

91 

455 

5 

15 

136 

816 

6 


One can see that for problems with more than 6 design variables, it will probably be more 
computationally efficient in an optimization algorithm to utilize a sequence of quadratic 
response surfaces than one 3rd or higher order response surface. When there are 6 or fewer 
design variables, 3rd or 4th order response surfaces may be used to advantage. 

In this report, the term "minimum point design" refers to a design that has just enough 
design points to allow the determination of coefficients of an approximating polynomial. 
The term "augmented minimum point design" is a minimum point design which contains 
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additional design points. Thus, augmented minimum point designs are over-determined 
designs. The studies that have been performed in this report indicate that augmented 
minimum point designs are competitive with, if not better than, central composite designs 
for developing a 2nd order response surface. A program DESIGNS [20] was developed for 
generating augmented minimum point designs for developing a 2nd order response surface. 
That program is described in Section 8.2. 

When there are 6 or fewer design variables, it may be computationally beneficial to use a 
3rd order or 4th order response surface. Thus, the program DESIGN4 [28] was developed 
to generate augmented minimum point designs for a 4th order response surface. The 
program DESIGN4 is discussed in Section 8.3. The program can also be used to develop 
a 3rd order response surface. The 3rd order minimum point design is a subset of the 4th 
order minimum point design. Thus the 4th order minimum point design will give an over- 
determined 3rd order approximation. Additional randomly selected design points can be 
added to the 4th order minimum point design to give the desired degree that the 3rd order 
approximation is to be over-determined. 

8.2 Augmented Minimum Point Designs for 2nd Order Approximations 

The basic building block for program DESIGNS is the star pattern of design points. Figure 
4 shows the star pattern for 3 design variables. This pattern of design points allows one to 
determine those coefficients of a 2nd order polynomial approximation associated with the 
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terms 


1, x p xf, i=M 


( 102 ) 


To be able to determine the coefficients associated with the terms 

xpc f i*j 


(103) 


one must supplement the star pattern with one additional design point in the Xj, Xj planes. 
Figure 30 shows the additional design point in the x^ Xj plane. Figure 6 shows the total 
minimum point design for 3 design variables. 


Studies of this report indicate that designs should be over-determined. Having a design that 
is 20-50% over-determined is a good compromise between keeping down the number of 
design points while still getting a good approximation. The program DESIGNS augments 
the minimum point design with a user selected number of random design points. 

8.2.1 Specifics of program DESIGNS 

A listing of the FORTRAN program DESIGNS is found in Appendix 1 and a copy of that 
program is found in file "designs.f on the floppy disk accompanying this report. The 
program should be compiled with a ¥11 compiler with the compiled program called "design". 
To run the program just enter "design" from the keyboard. The program prompts the user 
for 


80 



1. the number of design variables, 

2. the number of designs points to augment the minimum point design, and 

3. a seed parameter, IFLAG, which is used to generate the random numbers (IFLAG can 
be entered as any positive integer). 

The program then generates a design in local coordinates with the maximum range on each 
design variable of -1 to +1. The program then 

4. asks the user to enter an integer which specifies whether design point coordinates are 
to be also generated in global coordinates. If they are to be calculated in global 
coordinates, the program then 

5. prompts the user to enter the range of design variables in global coordinates. 

Results with commentary are written to file "design.res\ Design points without commentary 
are written to file "design.run". 


8.3 Augmented Minimum Point Design for 3rd and 4th Order Approximat ion 
A 3 k factorial design is used as the building block of this minimum point design. The 3 k 
factorial design provides information for calculating the coefficients associated with the 
terms 


1 , X., XpCjy Xi , Xj Xp 


2 2 
Xj Xj , 


J*l 


(104) 


Additional points are then added at -1 and 1 (in local coordinates) along the x ; axis. These 
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points together with the 3 k factorial design point give the star pattern which can be seen in 
Figure 31. With this arrangement of points, there are 5 design points along the x, axis which 
provides information for calculating the coefficient associated with the terms 


4 


( 105 ) 


Additional design points are then placed in each Xj, Xj plane which provides information for 
calculating the coefficient associated with the terms 


These points are also shown in Figure 31. 

8.3.1 Specifics of program DESIGN4 

A listing of the FORTRAN program DESIGN4 is found in Appendix 2 and a copy of that 
program is found in file "design4.f on the floppy disk accompanying this report. The 
program should be compiled with a F77 compiler with the compiled program called 
"design4". To run the program just enter "design4" from the keyboard. The program 
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prompts the user for needed information. Prompts and response are similar to those for the 
program DESIGNS. 

8.4 Conclusion 

A minimum point design is a design that has just enough design points to allow the 
determination of the coefficients of an approximating polynomial. An augmented minimum 
point design is a minimum point design which contains additional design points. Augmented 
mini mum point designs are competitive with, if not better than, central composite designs 
for developing a 2nd order response surface. Minimum point designs should be augmented 
with enough points that the approximation is 20-50% over-determined. A program 
DESIGNS was developed for generating augmented minimum point designs for developing 
a 2nd order response surface. 


When there are more than 6 design variables, 3rd or higher order approximations require 
so many design points that it is computationally better to perform a sequence of 2nd order 
approximations in an optimization process than one higher order approximation. When 
there are 6 or fewer design variables, a 2nd order approximation may often be satisfactory. 
However, for those cases where it is desirable to use a higher order approximation, program 
DESIGN4 was developed. That program generates designs which can be used to develop 
3rd or 4th order approximations. 
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9. Solution Algorithm 


9.1 Introduction 

In this investigation, the program NEWPSI was used to perform the studies involving 
polynomial approximations. That program can investigate under-determined, exactly- 
determined, or over-determined approximations of various orders. The version submitted 
with this report can handle up to 15 design variables as programmed. The order of 
polynomial it can handle is as follows: 

1. one design variable, up to a 20th order polynomial 

2. two design variables, up to a 5th order polynomial 

3. for 2-15 design variables, a second order polynomials. 

One can use up to 250 designs to train the approximation. In calculating v G , it can handle 
up to 2000 grid points. 

The program solves for the undetermined parameters associated with the approximation. 
It then evaluates the approximate function at the design points and calculates the error 
parameter, v. It then reads in the design points and function value on the test grid. The 
approximate function is evaluated at the grid points and the error parameter, v G , is then 
evaluated. 

9.2 Program Specifics 

A listing of the FORTRAN program NEWPSI is found in Appendix 3 and a copy of that 
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program is found in file "newpsi.f on the floppy disk accompanying this report. The 
program should be compiled with a F77 compiler and the compiled program called newpsi . 
To run the program just enter "newpsi" from the keyboard. Data is read from the file 
"newpsi.dat". Data can be in free format. The program asks for the following data: 

1. a value of the print code, ip; (If ip =4, great quantities of output are generated for 
program checkout. Normally the program is run with ip=0 for normal output). 

2. the number of design variable, nd; 

3. the order of the polynomial being considered, np; 

4. the number of design points in the design, m; 

5. the design and function value at the design points, x(i,j), y(i); 

6. the number of design points on the grid, ng; and 

7. the design and function value at the grid points, xx(i,j), yy(j)- 

Output is written to the screen and to file "newpsi.res". 
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10. Conclusion 


For a given order of approximation of a function, f, the quality of the approximation is 
affected by 

a. the number of levels of the design variables, 

b. the location of the design points, and 

c. the degree which the approximation is over-determined. 

For an nth order approximation, 

1. there must be n + 1 levels of the design variables; 

2. the design points must be located so that information is available for calculating all 
of the nth derivatives of f; 

3. the approximation should be, at least, 20-50% over-determined. 


For example, for a 2nd order approximation in 3 design variables, there must be at least 3 
levels of the design variables, design points must be located so that information is available 
for calculating 


jr 

dx’ dxpx’ 


i=l, 3; ;=1,3 


(107) 


A complete 2nd order polynomial approximation contains 10 undetermined coefficients. 
Thus, at least 10 design points are required to provide information for calculating these 
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coefficients. To have the approximation 30% over-determined, one would want to use 13 
design points. 

For second order approximations, when there are fewer than 6 design variables, central 
composite designs meet requirements 1-3. However, for 6 or more design variables, these 
designs contain too many design points. A minimum point design is one which contains just 
enough design points, meeting the derivative requirements of item 1 and 2 above, to exactly- 
determine the approximation. An augmented minimum point design is a minimum point 
design supplemented with additional design points. The program DESIGNS was developed 
to yield augmented minimum point designs for 2nd order approximations. The quality of 
approximations developed using designs from program DESIGNS was comparable to, if not 
better than, other standard designs such as the central composite designs. 

For more than 6 design variables, 3rd and 4th order approximations require so many design 
points to determine the coefficients in those approximations that it is more computationally 
efficient to develop a number of 2nd order approximations than one approximation of 3rd 
or higher order. For 6 or fewer design points, 2nd order approximations may be quite 
adequate. However, for those cases where one wishes to use a 3rd or 4th order 
approximation, the program DESIGN4 was develop. That program generates an augmented 
minimum point design for developing a 4th order approximation. 

Previous studies have shown that the quality of approximations using neural networks is 


87 


comparable to those using polynomial approximations when the number of undetermined 
parameters associated with the approximations is the same. Thus, neural networks trained 
with designs from DESIGNS or DESIGN4 should offer approximations of comparable 
quality to those obtained using polynomial approximations with the same number of 
associated undetermined parameters. 
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Figure 1. Artificial neural net 
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Figure 9. One dimensional example 
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Figure 10. 


D, A, E, and G optimality, first order approximation 
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Figure 11. D, A, E, and G optimality, second order approximation 
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D, A, E, and G optimality, third order approximation 




Parameter vG 

Fourth Order Polynomial Approximation 



Uniform — • — D Optimality A Optimality 

E Optimality — G Optimality 


Figure 13. D, A, E, and G optimality, fourth order approximation 
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Figure 14. S and Q optimality, first order approximation 
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Figure 15. S and Q optimality, second order approximation 
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Figure 16. S and Q optimality, third order approximation 
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Figure 17. S and Q optimality, fourth order approximation 
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Figure 18. Q optimality, 11 out of 13 points selected 
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Figure 19. Q optimality, 7 out of 13 points selected 
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Figure 20. Q optimality, 5 out of 13 points selected 
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Figure 21. Random points, first order approximation 
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Figure 22. Random points, second order approximation 
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Figure 23. Random points, third order approximation 
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Figure 24. Random points, fourth order approximation 
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Figure 25. Significance testing, Example 1, 2nd order 
approximation 
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Significance testing. Example 1, 3rd order 
approximation 
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Appendix 1 

Program DESIGNS 

PROGRAM DESIGNS 
C 

C PROGRAM TO GENERATE DESIGNS FOR 2ND ORDER POLYNOMIAL 

-C PROGRAM DIMENSIONED FOR UP TO 20 VARIABLES 

C RESULTS TO SCREEN AND TO FILE designs. res 

C DESIGN IN GLOBAL COORDINATES TO FILE designs. run 

_C 

C DEFINITIONS 

C N = NUMBER OF DESIGN VARIABLES 

CM = NUMBER OF RANDOM DESIGNS POINTS 

“C 

DIMENSION X (2000 ,20) 

DIMENSION XBB (20) ,XBE(20) ,A(20) ,B(20) 

1 FORMAT (I5,6F10.6) 

2 FORMAT (' PROGRAM GENERATES DESIGNS FOR FITTING 2ND ORDER', 

X' POLYNOMIAL') 

3 FORMAT (' ENTER NUMBER OF DESIGN VARIABLES') 

4 FORMAT (' NUMBER OF DESIGN VARIABLES = N =' , 13) 

11 FORMAT (6F10. 6) 

OPEN (UNIT=7 , FILE= ' designs . res ' ) 

OPEN (UNIT=8 , FILE= ' designs . run ' ) 

WRITE (6,2) 

WRITE (6,3) 

READ (5, *)N 
WRITE ( 6 , 4 ) N 
SET UP TERMS 
NP1=N+1 
NM1=N-1 

M=(N*N+3*N+2) /2 
MP1=M+1 

ZERO DESIGN MATRIX 
D0100I=1 , M 
D0100J=1 , N 
100 X(I,J)=0. 

11=0 



~c 

C GENERATE THE FIRST N+l POINTS FOR FITTING A LINEAR FUNCTION 

C THE FIRST POINT IS WHEN ALL X'S ZERO, ALREADY DONE 

-C GENERATE NEXT N POINTS 

DO101I=l,N 
11=1+1 

_ 101 X(II,I)=1. 

C 



C GENERATE NEXT N POINTS 

~C THE 2N+1 POINTS THUS GENERATED WILL ALLOW ADDING SQUARED TERMS 
C 

D0102I=1 , N 
II=I+N+1 
102 X(II,I)=-1. 

C 

_c 

c 

C GENERATE NEXT N(N-l)/2 POINTS 

C THE (N*N+3*N+2 ) /2 POINTS THUS GENERATED WILL ALOW ADDING CROSS 

— c PRODUCT TERMS. WE WILL THEN HAVE COMPLETE 2ND ORDER POLYNOMIAL 

C APPROXIMATION 

C 



on ooooo non non 


ILAST=2*N+1 

IDO=N-l 

J=1 

JJ=2 

103 CONTINUE 
DO104I=l , IDO 
II=I+ILAST 
X(II, J)=l. 

X(II, JJ)=1. 

JJ=JJ+1 

104 CONTINUE 
ILAST=ILAST+IDO 
IDO=IDO-l 
J=J+1 

JJ=J+1 

IF ( J. LE. NM1) GOTO 10 3 

IF WE GOT HERE WE HAVE DEVELOPED THE MINIMUM POINT DESIGN 

WRITE ( 6 , * ) ' WE HAVE GENERATED ',11, ' POINTS IN THE MIN PT DESIGN' 
WRITE (7 , *) ' WE HAVE GENERATED ',11,' POINTS IN THE MIN PT DESIGN' 
WRITE ( 6 , * ) ' DESIGN POINTS WRITTEN TO FILE designs. res' 

DEVELOP DESIGN POINTS TO AUGMENT THE MINIMUM POINT DESIGN 
READ IN THE NUMBER OF RANDOM DESIGN POINTS TO BE DEVELOPED 
WRITE (6,*)' ENTER THE NUMBER OF RANDOM GENERATED DESIGN PTS', 

X' DESIRED=M' 

READ ( 5 , * ) M 

WRITE (6 , *) ' NUMBER OF RANDOM DESIGN POINTS M=' ,M 
WRITE (7,*)' NUMBER OF RANDOM DESIGN POINTS M=',M 

WRITE ( 6 , * ) ' I FLAG IS ANY POSITIVE INTEGER USED TO START RANDOM', 
X' PROCESS' 

WRITE (6,*)' ENTER I FLAG' 

READ ( 5 , * ) I FLAG 
WRITE ( 6 , * ) ' IFLAG= ' , I FLAG 
WRITE ( 7 , * ) ' IFLAG= ' , I FLAG 
D0850I=1,M 
11=11+1 
D0851J=1 , N 
IFLAG=IFLAG+1 
XDUM=RAND ( I FLAG ) 

X(II, J)=2.*XDUM-1. 

851 CONTINUE 
850 CONTINUE 

IF WE GOT HERE WE HAVE FINISHED GENERATING THE RANDOM DESIGN PTS 
WRITE ( 6 , * ) ' RANDOM DESIGN POINTS WRITTEN TO FILE designs. res' 

PRINT OUT THE MINIMUM POINT MATRIX IN LOCAL COORDINATES 

WRITE ( 7 , * ) ' DESIGN MATRIX IN LOCAL COORDINATES' 

ITOTAL=II 

D0700I=1 , ITOTAL 

WRITE (7 ,1)1, (X(I,J) , J=1 , N) 

700 CONTINUE 

SEE IF WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 

WRITE (6 , * ) ' ITEST=1 IF DESIGN POINTS ARE TO BE IN GLOBAL', 

X' COORDINATES' 

WRITE (6,*)' OTHERWISE, ITEST=0' 



non 


WRITE (6,*)' ENTER ITEST' 

READ ( 5 , * ) ITEST 
IF (ITEST. NE. 1) GOTO860 

-C IF WE GOT HERE WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 

WRITE ( 6 , * ) ' ENTER LOWER AND UPPER RANGE ON EACH DESIGN VARIABLE ' 
WRITE (6 , *) ' i.e. ENTER XBB(I) TO XBE(I)' 

D0861I=1 , N 

READ ( 5 , * ) XBB ( I ) , XBE ( I ) 

WRITE (6,*) 7 I, XBB (I) , XBE (I) =' , I , XBB (I) ,XBE(I) 

WRITE (7 , *) ' I, XBB (I) , XBE ( I ) = ' , I , XBB ( I ) ,XBE(I) 

- 861 CONTINUE 
G0T0862 

860 CONTINUE 

IF WE GOT HERE LOWER BOUND VARIABLE IN GLOBAL COORDINATES IS -1 
IF WE GOT HERE UPPER BOUND VARIABLE IN GLOBAL COORDINATES IS 1 
D0863I=1 , N 
XBB ( I ) =-l . 

XBE (I ) =1 . 

863 CONTINUE 

- 862 CONTINUE 

WRITE (7, *) ' I, XBB (I) ,XBE(I) ,A(I) ,B(I) ' 

D01301I=1 , N 

A ( I ) = ( XBE ( I ) -XBB ( I ) ) / 2 . 

B ( I ) = ( XBE ( I ) +XBB ( I ) ) / 2 . 

WRITE (7, *) I, XBB (I) ,XBE(I) ,A(I) ,B(I) 

1301 CONTINUE 

- D012 02 1=1 , ITOTAL 
DO1202 J=1 , N 

1202 X(I, J)=A(J) *X(I, J) +B(J) 

- WRITE ( 6 , * ) ' DESIGN IN GLOBAL COORDINATES WRITEN TO designs. res' 
WRITE(6, *) ' DESIGN IN GLOBAL COORDINATES WRITEN TO designs. run' 
WRITE (7,*)' DESIGN IN GLOBAL COORDINATES' 

__ WRITE ( 8, *) ITOTAL 

D097 01=1, ITOTAL 
WRITE (7,1)1, (X(I,J) , J=1 , N) 

WRITE (8,11) (X(I,J) , J=1 , N) 

- 970 CONTINUE 

STOP 

END 
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Appendix 2 

Program DESI6N4 

PROGRAM DESIGN4 
C 

C PROGRAM TO GENERATE DESIGNS FOR 4TH ORDER POLYNOMIAL 

-C PROGRAM DIMENSIONED FOR UP TO 6 VARIABLES 

C RESULTS TO SCREEN AND TO FILE design4.res 

C DESIGN IN GLOBAL COORDINATES TO FILE design4.run 

_C 

C DEFINITIONS 

C N = NUMBER OF DESIGN VARIABLES 

CM = NUMBER OF RANDOM DESIGNS POINTS 

~C 

DIMENSION X (2000, 6) 

DIMENSION XBB(IO) ,XBE(10) ,A(10) ,B(10) 

1 FORMAT (15 , 6F10. 6) 

2 FORMAT ( ' PROGRAM GENERATES DESIGNS FOR FITTING 4TH ORDER', 

X' POLYNOMIAL') 

3 FORMAT (' ENTER NUMBER OF DESIGN VARIABLES') 

4 FORMAT ( ' NUMBER OF DESIGN VARIABLES = N =', 13) 

11 FORMAT (6F10. 6) 

OPEN (UNIT=7 , FILE= ' design4 . res ' ) 

OPEN (UNIT=8 , FILE= ' design4 . run ' ) 

WRITE (6,2) 

WRITE (6,3) 

READ ( 5 , * ) N 
WRITE ( 6 , 4 ) N 
IF(N.EQ. 6)GOTO601 
IF (N. EQ. 5) GOTO501 
IF (N . EQ . 4 ) G0T04 0 1 
IF (N . EQ. 3 ) G0T03 01 
IF (N. EQ. 2 ) G0T02 01 
IF (N . EQ. 1) GOTOlOl 

WRITE ( 6 , * ) ' PROGRAM CAN NOT DO MORE THAN 6 DESIGN VARIABLES ' 

WRITE ( 7 , * ) ' PROGRAM CAN NOT DO MORE THAN 6 DESIGN VARIABLES' 

STOP 

DEVELOP 3 FACTORIAL DESIGN TO GET 4 DESIGN VARIABLE PRODUCT TERMS 
101 CONTINUE 
11=0 

D0100I1=1, 101, 50 
11=11+1 

X (II, l)=FLOAT( 11-51) / 100. 

100 CONTINUE 
GOTO701 
201 CONTINUE 
11=0 

D0200I1=1, 101,50 
D0200I2=1 ,101,50 
11=11+1 

X (II, l)=FLOAT ( 11-51) / 100. 

X (II , 2 ) =FLOAT (12-51) / 100 . 

_ 200 CONTINUE 

GOTO701 
301 CONTINUE 
11=0 

D0300I1=1, 101,50 
D03 0012=1, 101,50 
D0300I3=1 ,101,50 
11=11+1 

X ( II , 1) =FLOAT (11-51) / 100. 

X ( I I , 2 ) =FLOAT (12-51) /100. 



non 


X ( II , 3 ) =FLOAT (13-51) /100. 

300 CONTINUE 
GOTO701 

- 401 CONTINUE 

11=0 

D0400I1=1, 101, 50 
_ D0400I2=1, 101, 50 

D0400I3=1, 101, 50 
DO400I4=l, 101,50 
11=11+1 

X (II, l)=FLOAT( 11-51 )/ 100. 

X (II , 2 ) =FLOAT (12-51) /100. 

X(II, 3 ) =FLOAT ( 13-51 ) /100. 

- X (II , 4 ) =FLOAT (14-51) /100. 

400 CONTINUE 

GOTO701 
501 CONTINUE 
11=0 

DO500Il=l, 101, 50 
D0500I2=1, 101,50 
~ DO500I3=l , 101,50 

D0500I4=1 , 101,50 
D0500I5=1, 101,50 
_ 11=11+1 

X (II , 1) =FLOAT ( 11-51) /100. 

X (II , 2 ) =FLOAT (12-51) /100. 

X (II , 3 ) =FLOAT ( 13-51) /100. 

X (II , 4 ) =FLOAT (14-51) /100. 

X (II , 5) =FLOAT (15-51) /100. 

500 CONTINUE 

- GOTO701 
C 

601 CONTINUE 
_ 11=0 

D0600I1=1 ,101,50 
D0600I2=1, 101, 50 
D0600I3=1 ,101,50 
_ D0600I4=1 ,101,50 

DO600I5=l ,101,50 
D0600I6=1, 101, 50 

- 11=11+1 

X ( II , 1) =FLOAT ( 11-51 ) /100. 

X (II , 2 ) =FLOAT ( 12-51 ) /100. 

_ X (II , 3 ) =FLOAT ( 13-51 ) /100. 

X (II , 4 ) =FLOAT ( 14-51 ) /100. 
X(II,5)=FLOAT(I5-51) /100. 
X(II,6)=FLOAT(I6-51) /100. 

~ 600 CONTINUE 

GOTO701 

701 CONTINUE 

ENTER REST OF POINTS IN THE STAR FORMATION 

D0702I=1 , N 
11=11+1 
DO703 J=1 , N 
703 X(II,J)=0. 

X(II,I)=1. 

702 CONTINUE 
DO704I=l , N 
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11=11+1 
D0705J=1,N 
705 X(II,J)=0. 

- X(II,I)=-1. 

704 CONTINUE 

C 

_C ENTER TERMS TO CALCULATE COEFFICIENT ASSOCIATED WITH THE TERM 

C X(I) **3*X ( J) 

C 

NM1=N-1 
~ IDO=N-l 

J=1 
JJ=2 

- 803 CONTINUE 

D0804I=1, IDO 
11=11+1 
X(II,J)=1. 

X ( II , JJ) = . 5 
11=11+1 
X(II, J)=.5 

- X(II,JJ)=1. 

JJ=JJ+1 

804 CONTINUE 
IDO=IDO-l 
J=J+1 
JJ=J+1 

IF ( J. LE. NM1) GOTO803 

IF WE GOT HERE WE HAVE DEVELOPED THE MINIMUM POINT DESIGN 

WRITE ( 6 , * ) ' WE HAVE GENERATED ' ,11,' POINTS IN THE MIN PT DESIGN 
WRITE ( 7 , * ) ' WE HAVE GENERATED ',11,' POINTS IN THE MIN PT DESIGN 
WRITE (6,*)' DESIGN POINTS WRITTEN TO FILE design4.res' 

DEVELOP DESIGN POINTS TO AUGMENT THE MINIMUM POINT DESIGN 
READ IN THE NUMBER OF RANDOM DESIGN POINTS TO BE DEVELOPED 
WRITE (6,*)' ENTER THE NUMBER OF RANDOM GENERATED DESIGN PTS', 

X' DESIRED=M' 

READ ( 5 , * ) M 

WRITE (6,*)' NUMBER OF RANDOM DESIGN POINTS M=' ,M 
WRITE (7 , *) ' NUMBER OF RANDOM DESIGN POINTS M=' ,M 

WRITE (6,*)' I FLAG IS ANY POSITIVE INTEGER USED TO START RANDOM', 
X' PROCESS' 

WRITE (6,*)' ENTER I FLAG' 

READ ( 5 , * ) I FLAG 
WRITE ( 6 , * ) ' IFLAG= ' , I FLAG 
WRITE (7 , *) ' IFLAG= ' , I FLAG 
DO850I=l,M 
11=11+1 
D0851J=1 , N 
IFLAG=IFLAG+1 
XDUM=RAND ( I FLAG ) 

X(II, J)=2.*XDUM-1. 

851 CONTINUE 
850 CONTINUE 

IF WE GOT HERE WE HAVE FINISHED GENERATING THE RANDOM DESIGN PTS 
WRITE (6,*)' RANDOM DESIGN POINTS WRITTEN TO FILE design4.res' 

PRINT OUT THE MINIMUM POINT MATRIX IN LOCAL COORDINATES 
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c 


WRITE (7 , *) ' DESIGN MATRIX IN LOCAL COORDINATES' 

ITOTAL=II 

D07 001=1, ITOTAL 

WRITE (7 ,1)1, (X(I,J) , J=1 , N) 

700 CONTINUE 

SEE IF WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 

WRITE ( 6 , * ) ' ITEST=1 IF DESIGN POINTS ARE TO BE IN GLOBAL', 

X' COORDINATES' 

WRITE ( 6 , * ) ' OTHERWISE, ITEST=0' 

WRITE (6,*)' ENTER ITEST' 

READ (5, *) ITEST 
IF (ITEST. NE . 1) GOTO860 

IF WE GOT HERE WE ARE TO GENERATE DESIGNS IN GLOBAL COORDINATES 
WRITE (6,*)' ENTER LOWER AND UPPER RANGE ON EACH DESIGN VARIABLE' 
WRITE ( 6 , * ) ' i.e. ENTER XBB(I) TO XBE(I)' 

D0861I=1 , N 

READ ( 5 , * ) XBB ( I ) ,XBE(I) 

WRITE (6,*) ' I, XBB (I) , XBE ( I ) = ' , I , XBB ( I ) ,XBE(I) 

WRITE (7,*) ' I, XBB (I) , XBE ( I ) = ' , I , XBB (I ) ,XBE(I) 

861 CONTINUE 
G0T0862 

860 CONTINUE 

IF WE GOT HERE LOWER BOUND VARIABLE IN GLOBAL COORDINATES IS -1 
IF WE GOT HERE UPPER BOUND VARIABLE IN GLOBAL COORDINATES IS 1 
D0863I=1 , N 
XBB (I ) =-l . 

XBE (I ) =1 . 

863 CONTINUE 

862 CONTINUE 

WRITE ( 7 , * ) ' I , XBB (I) ,XBE(I) ,A(I) ,B(I) ' 

DO1301I=l , N 

A (I) = (XBE (I) -XBB (I) ) /2. 

B (I) = (XBE (I) +XBB (I) ) / 2 . 

WRITE ( 7 , * ) I , XBB ( I ) ,XBE(I) ,A(I) ,B(I) 

1301 CONTINUE 

D01202I=1, ITOTAL 
D01202 J=1 , N 

1202 X(I, J)=A(J) *X(I, J)+B(J) 

WRITE (6 , *) ' DESIGN IN GLOBAL COORDINATES WRITEN TO design4.res' 
WRITE (6 , *) ' DESIGN IN GLOBAL COORDINATES WRITEN TO design4.run' 
WRITE (7,*)' DESIGN IN GLOBAL COORDINATES' 

WRITE (8 , *) ITOTAL 

D0970I=1, ITOTAL 

WRITE (7, 1) I, (X(I, J) , J=1,N) 

WRITE(8 , 11) (X(I,J) ,J=1,N) 

970 CONTINUE 
STOP 
END 
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Appendix 3 
Program NEWPSI 


PROGRAM newpsi 

— c 

c ************************************************************* 

c ************************************************************* 

— c 

c the program develops a polynomial approximation which 

c may be either under, exactly, or over determined 

c it can handle up to 15 design variables as programmed. 

c The order of polynomial it can handle is as follows: 

c 1. one one design variable, up to a 20th order polynomial 

c 2. two design variables, up to 5th order polynomial 

~ c 3- for 2-15 design variables, a 2nd order polynomial 

c One can use up to 250 designs to train the approximation, 

c It can handle up to 2000 grid points 

_c 

c ************************************************************ 

c ************************************************************ 

c 

~ IMPLICIT REAL *8 (A-H,0-Z) 

dimension x (250, 15) ,y (250) , a (250, 136) 
dimension yhat(250) 

— dimension b(136) 

dimension xx(2000, 15) ,yy(2000) ,abig(2000, 136) 
dimension yyhat(2000) 

_ 1 FORMAT (9F8. 4) 

2 FORMAT ( 3 F12 . 6) 

3 FORMAT (F10. 6, 1H, ,F10.6,1H, ,F10.6,1H, ,F10.6,1H, ,F10.6,1H, ,F10.6, 
X1H, , F10 . 6) 

“ OPEN (UNIT=5,FILE= / newpsi.dat') 

OPEN (UNIT=7 , FILE= ' newpsi . res' ) 

OPEN (UNIT=8 , FILE= ' newpsi . plot ' ) 

-C 

c *************************************************************** 

read in data 
c 

c read in the print code 

read (5, *) ip 
'c 

c enter number of design variables, nd 

read (5, *) nd 

c enter THE DEGREE OF POLYNOMIAL TO BE CONSIDERED, np 

READ ( 5 , * ) np 

c ENTER NUMBER OF DESIGNS FOR PROBLEM, M 

READ ( 5 , * ) M 

write (6,*)' print code ip=',ip 

write(6,*)' number of design variables, nd=',nd 



noon 


write(6,*)' degree of polynomial being considered-np-' , np 
write (6,*) ' number of designs m=',m 
write (7,*)' print code ip=',ip 

— write (7,*) ' number of design variables, nd— ,nd 
write(7,*)' degree of polynomial being considered=np=' , np 
write (7,*) ' number of designs m=',m 

_C 

c read in designs and set up matrix a 

C 

write (7 , *) ' x(i, j) ,y(i) ' 

- D0101I=1,M 

read (5, *) (x(i, j) , j=l,nd) ,y(i) 
write(7,*) (x(i, j) , j=l,nd) ,y(i) 

_ 101 continue 

c set up the coefficient matrix, a, in the matrix equation 

_c y=a x 

c 

call geta (ip,m, nd, np, n, x, a) 

— C SEE WHETHER SYSTEM IS UNDER, EXACTLY , OR OVER DETERMINED 

C 

IF (M.GE.N) GOT0400 

_C IF WE GOT HERE WE ARE UNDER-DETERMINED 

WRITE ( 6 , * ) ' SYSTEM IS UNDER-DETERMINED ' 

WRITE (7 , *) ' SYSTEM IS UNDER-DETERMINED ' 

CALL PSI (ip,M,N,A, Y,B) 

GOTO402 

400 CONTINUE 

IF (M . GT . N ) G0T04 0 1 

— C IF WE GOT HERE WE ARE EXACTLY DETERMINED 

WRITE (6, *) ' SYSTEM IS EXACTLY DETERMINED' 

WRITE (7,*)' SYSTEM IS EXACTLY DETERMINED' 

__ CALL EXACT ( ip , M , A , Y , B ) 

GOTO402 

401 CONTINUE 

C IF WE GOT HERE WE ARE OVER-DETERMINED 

— WRITE(6,*)' SYSTEM IS OVER-DETERMINED' 

WRITE (7 , * ) ' SYSTEM IS OVER-DETERMINED' 

CALL OVER (ip,M,N,A,Y,B) 

- 402 CONTINUE 


EVALUATE APPROXIMATION AT DESIGNS 

WRITE ( 6 , *) ' MATRIX OF COEFFICIENTS, B(I)' 

WRITE (7 , *) ' MATRIX OF COEFFICIENTS, B(I)' 

WRITE ( 6 , * ) (B (I ) ,1=1, N) 

WRITE ( 7 , * ) (B (I) ,1=1, N) 

WRITE (7 , *) ' MATRICES Y(I) AND YHAT ( I ) ' 

_c 

c recalculate matrix a 

call geta(ip,m,nd,np,n,x,a) 

~c calculate approximation at designs and print results 

c 

write (7 , *) ' y (i) ,yhat (i) ' 

D0102I=1,M 
YHAT (I) =0. 

DO103 J=1 , N 



103 

- 102 
C 
c 


- 601 


~~ 603 


c 
— c 

c 

602 

c 

c 

c 

c 


c 

c 

c 

c 

c 


— c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


yhat ( i ) =yhat (i)+a(i,j)*b(j) 

CONTINUE 

WRITE (7,*)Y(I) , YHAT ( I ) 

CONTINUE 

evaluate function at grid 
read(5,*)ng 

write(6,*)' number of designs on grid = ngn',ng 
write (7,*)' number of designs on grid = ngn',ng 
write (7,*)' xx (i, j) ,yy(i) ' 

DO601I=l , ng 

read (5,*) (xx(i, j) , j=l,nd) ,yy(i) 
write (7 , *) (xx(i, j) , j=l,nd) ,yy(i) 
continue 

call getabg ( ip , ng , nd , np , n , xx , abig ) 
write (7 , *) ' yy(i) , yyhat (i) at grid' 

DO602I=l , ng 
YYHAT (I) =0 . 

D0603 J=1 , N 

yyhat ( i ) =yyhat ( i ) +abig ( i , j ) *b ( j ) 

CONTINUE 

WRITE ( 7 , * ) Y Y ( I ) , YYHAT (I) 
write the plot file 

write (8,*) (xx(i, j) , j=l,nd) , yyhat (i) 

CONTINUE 

calculate statistical terms 

call stat it (m , y , yhat , ng , yy , yyhat ) 

STOP 

END 

subroutine geta ( ip , m , nd , np , n , x , a ) 

************************************************************* 

************************************************************* 

This subroutine generates the matrix a where the matrix 
equation is y= a b. Here y are the training functions, 
b are undetermined coefficients. The algorithm is programmed 
to handle 

1. any level of approximation for one design variable 

2. up to 5th order polynomial in two design variables 

3. quadratic approximation in more than two design variabaales 

************************************************************ 

************************************************************ 

IMPLICIT REAL* 8 (A-H,0-Z) 
dimension x(250, 15) ,a(250, 136) 

do for each design 

do300i=l,m 

******************************************************** 
if nd is not equal to 1 go to 400 



if (nd.ne. I)goto400 


c 

c 

c 

c 

c 

c 

c 


201 


c 

c 

c 

400 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


******************************************************** 

******************************************************** 

here we have nd=l, i.e. one design variable 
we will develp a's for all np's 

a (i/ 1) =1. 
j=l 

do201k=l , np 
j=j+l 

a(i, j)=x(i, 1) **k 

continue 

n=np+l 

goto301 


************************************ 

continue 

if nd is not equal to 2 go to 500 
if (nd.ne. 2) goto500 

************************************************************** 

************************************************************** 

if we got here we have 2 design variables 

xl=x(i, 1) 
x2=x(i,2) 

********************* 

add the constant and linear terms 

a(i,l)=l. 

a(i,2)=xl 

a(i,3)=x2 

n=3 

if (np. It. 2) goto301 
******************** 

add the 2nd order terms 

a(i,4)=xl**2 

a(i,5)=xl*x2 

a(i,6)=x2**2 

n=6 

if (np. It . 3 ) goto301 
******************* 

add the cubic terms 

a(i,7)=xl**3 
a(i,8)=xl**2*x2 
a (i, 9) =xl*x2**2 
a(i, 10) =x2**3 



n=10 

if (np. It . 4) goto301 
c 

c ******************* 

c 

c add the 4th order terms 

c 

a (i , 11) =xl**4 
a(i, 12)=xl**3*x2 
a (i , 13) =xl**2*x2**2 
a(i, 14)=xl*x2**3 
a(i, 15)=x2**4 
n=15 

if (np. It . 5) goto301 
c 

Q ******************* 

C 

c add the 5th order terms 

c 

a (i, 16) =xl**5 
a(i, 17)=xl**4*x2 
a (i, 18) =xl**3*x2**2 
a (i, 19) =xl**2*x2**3 
a (i , 20) =xl*x2**4 
a(i,21)=x2**5 
n=21 

if (np. It. 6 ) goto3 01 
c 

c ****************** 

c 

c algorithm not programed for polynomials of order larger than 5 

c 

write(6,*) / for two design variables, algorithm not programed for' 
write (6,*)' polynomials of order larger than 5' 

write(7,*)' for two design variables, algorithm not programed for' 
write(7,*)' polynomials of order larger than 5' 
stop 
c 

c ****************************************************** 

c ****************************************************** 

c 

500 continue 


c 

c 

c 

c 

c 

c 

c 


501 


c 

c 

c 


if we got here number of design variables >2 
******************** 


enter constant and linear terms 

a(i,l)=l. 

j=l 

do501k=l,nd 

j=j+l 

a(i, j)=x(i,k) 

continue 

n=j 

if (np. It. 2) goto301 
******************** 



c enter the quadratic terms 

c 

do502k=l,nd 
do502L=k, nd 
j=j+i 

a(i, j)=x(i,k) *x(i,L) 

502 continue 
n=j 

if (np. It. 3) goto301 
c 

c ******************** 

c algorithm not programmed for more than quadratic approximation 

c when number of design variables >2 

write (6, *) ' algorithm not programmed for more than quadratic' 
write(6,*)' approximation when number of design variables >2^ 
write(7,*)' algorithm not programmed for more than quadratic 
write (7,*)' approximation when number of design variables >2 

stop 


— c 
c 
c 
_c 
c 
c 




********* 

********* 


** 

* 


print out some results 

301 continue 

if (ip. It. 4) goto302 

write (6,*) ' a(i,j) (a(i,j) , j=l,n) 
write (6 ,*) ' ' 

write (7,*) ' a(i,j)',(a(i,j),j=l,n) 
write (7 , *) ' ' 

302 continue 


c 

c 

c 




300 continue 

wr ite(6,*)' number of undetermined coef=n- , n 
write(7,*)' number of undetermined coef=n=',n 


return 

end 

subroutine getabg ( ip , m, nd , np , n , x , a) 


c 
" c 
c 
c 

“C 
C 
C 
— C 

c 


_c 
— c 
c 
c 
— c 


************** 

************** 


*********************************************** 

*********************************************** 


This subroutine generates the matrix a where the matrix 
equation is y= a b. Here y are the training functions, 
b are undetermined coefficients. The algorithm is programmed 

to handle , . _ 

1. any level of approximation for one design variable 

2. up to 5th order polynomial in two design variables 

3. quadratic approximation in more than two design vanabaales 


********* 

********* 


*************************************************** 

*************************************************** 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A ( 2000 , 136) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


201 


c 

c 

c 

400 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


DIMENSION X ( 2000 , 15) 
do for each design 
do300i=l,m 

******************************************************** 

if nd is not equal to 1 go to 400 
if (nd.ne. 1) goto400 

******************************************************** 

******************************************************** 

here we have nd=l, i.e. one design variable 
we will develp a's for all np's 

a ( i , 1) =1 . 
j=l 

do201k=l , np 

j=j+l 

a(i, j)=x(i,l)**k 

continue 

n=np+l 

goto301 


************************************ 

continue 

if nd is not equal to 2 go to 500 
if (nd.ne. 2) goto500 

************************************************************** 

************************************************************** 

if we got here we have 2 design variables 

xl=x(i, 1) 
x2=x(i, 2) 

********************* 

add the constant and linear terms 

a (i, 1) =1 • 
a(i,2)=xl 
a(i,3)=x2 
n=3 

if (np. It. 2) goto301 

******************** 

add the 2nd order terms 

a(i,4)=xl**2 
a(i,5)=xl*x2 
a (i , 6) =x2**2 
n=6 

if (np. It . 3 ) goto301 



******************* 


c 
c 
c 

~c add the cubic terms 

c 

a (i , 7) =xl**3 
_ a (i, 8) =xl**2*x2 

a (i , 9) =xl*x2**2 
a(i, 10)=x2**3 
n=10 

if (np. It . 4) goto301 
c 

Q ******************* 

— C 

c add the 4th order terms 

c 

_ a(i f ll)=xl**4 

a (i, 12) =xl**3*x2 
a(i, 13) =xl**2*x2**2 
a (i , 14) =xl*x2**3 
“ a(i, 15)=x2**4 

n=15 

if (np. It. 5)goto301 
— c 

Q ******************* 

c 

_c add the 5th order terms 

c 

a (i, 16) =xl**5 
a (i, 17) =xl**4*x2 

- a(i, 18) =xl**3*x2**2 
a (i, 19) =xl**2*x2**3 
a (i, 20)=xl*x2**4 

_ a (i, 21)=x2**5 

n=21 

if ( np . It . 6 ) goto3 0 1 

_c 

Q ****************** 

c 

c algorithm not programed for polynomials of order larger than 5 

~c 

write(6,*)' for two design variables, algorithm not programed for' 

write (6,*)' polynomials of order larger than 5' 

_ write(7,*)' for two design variables, algorithm not programed for' 

write(7,*)' polynomials of order larger than 5' 

stop 
c 

~c ****************************************************** 

c ****************************************************** 

c 

— 500 continue 
c 

c if we got here number of design variables >2 

_c 

Q ******************** 

c 

c enter constant and linear terms 

— c 


a(i,l)=l. 

j=l 



do501k=l,nd 

j=j+l 

a(i, j)=x(i,k) 

501 continue 
n=j 

if (np. It . 2) goto301 
c 

Q ******************** 

c 

c enter the quadratic terms 

c 

do502k=l , nd 
do502L=k, nd 
j=j+l 

a(i, j)=x(i,k) *x (i , L) 

502 continue 

n=j 

if (np. It . 3 ) goto301 
c 

c ******************** 

c 

c algorithm not programmed for more than quadratic approximation 

c when number of design variables >2 

c 

write(6,*)' algorithm not programmed for more than quadratic' 
write(6,*)' approximation when number of design variables >2' 
write(7,*)' algorithm not programmed for more than quadratic' 
write(7,*)' approximation when number of design variables >2' 
stop 
c 

c ***************************************************************** 

c **************************************************************** 

c 

c print out some results 

c 

301 continue 

if (ip. It. 4) goto302 

write (6,*) ' a(i,j)',(a(i,j) , j=l,n) 
write ( 6 , * ) ' ' 

write (7 , *) ' a (i, j ) ', (a (i , j ) , j=l , n) 
write (7,*)' ' 

302 continue 
c 

c ******************************************************************* 

c 

300 continue 

write(6,*)' number of undetermined coef=n=',n 
write(7,*)' number of undetermined coef=n=',n 
c 

return 

end 

SUBROUTINE PSI (IP, M, N , DUMA, Y , XX) 

IMPLICIT REAL* 8 (A-H,0-Z) 

C 

DIMENSION DUMa (250,136) 

DIMENSION A(21, 21) ,B(21,21) ,D(21,21) ,DI(21,21) ,BPI(21,21) 

DIMENSION C(21,21) , FI (21,21) , CPI (21, 21) ,H(21,21),HI(21,21) 

DIMENSION API (21,21) 

DIMENSION F ( 2 1 , 21) 

DIMENSION IPIVOT (21), IWK (21,2) 



DIMENSION y(250) 

DIMENSION XX (21) 

—q THIS SUBROUTINE CALCULATES PSEUDO INVERSE OF MATRIX A 

c M ROW DIMENSION OF A LESS THAN N 

C N = COLUMN DIMENSION OF A 

_C 

C COPY DUMA TO A 

C 

D090I=1,M 
~ DO90J=l , N 

90 A ( I , J) =DUMA ( I , J) 

C 

— c 

c PRINT MATRIX A 

if ( ip. It . 4 ) gotoSO 
WRITE ( 6 , * ) ' MATRIX A' 

WRITE (7 , * ) ' MATRIX A' 

CALL WRITIT (M, N , A) 

50 continue 

-~c 

C SET UP MATRIX B 

C 

_ doiooi=i,m 

D0100J=1,M 

100 B(I, J) =A(I , J) 

if ( ip . It . 4 ) goto5 1 
" WRITE (6, *) ' MATRIX B' 

WRITE (7 , * ) ' MATRIX B' 

CALL WRITIT (M,M,B) 

— 51 continue 

C 

C GET D= B TRAN * B 

_C 

doioii=i,m 
doioij=i,m 
D (I , J) =0 . 

D0101K=1,M 

101 D (I , J) =D (I , J) +B (K, I) *B (K, J) 
if (ip. It. 4) goto52 

WRITE ( 6 , * ) ' MATRIX D' 

WRITE (7 , *) ' MATRIX D' 

CALL WRITIT (M,M,D) 

_ 52 continue 

C 

C GET INVERSE OF D=DI 

MAX=21 
“ MDUM=0 

CALL°MATINV (MAX ,M,D, MDUM , DI , IOP , DETERM , ISCALE , IPIVOT , IWK) 
- WRITE ( 6 , * ) ' DETERM =/ , DETERM , ' ISCALE=' , ISCALE 

WRITE ( 7 » * ) ' DETERM= ' , DETERM , ' I SCALE= / , ISCALE 

D0300I=1,M 
_ D0300J=1 , M 

300 DI (I , J) =D (I / J) 

if (ip. It. 4) goto53 
WRITE ( 6 , * ) ' MATRIX DI' 

WRITE ( 7 , * ) ' MATRIX DI' 

CALL WRITIT (M,M,DI) 

53 continue 


c 

c GET PSEUDO INVERSE OF B = BPI = DI * B TRANS 

C 

- DO102I=l , M 
D0102 JQ=1 , M 
BPI (I, JQ)=0. 

_ D0102 J=1 , M 

102 BPI (I , JQ) =BPI (I , JQ) +DI (I,J)*B(JQ,J) 
if (ip. lt.4)goto54 

WRITE ( 6 , * ) ' MATRIX BPI' 

“ WRITE (7 ,*) ' MATRIX BPI ' 

CALL WRITIT (M, M, BPI ) 

54 continue 

-c 

C SET UP MATRIX C = A 

C 

_ D0103I=1 , M 

D0103 J=1 , N 

103 C(I, J)=A(I, J) 

if (ip. It . 4) goto55 
~ WRITE (6 ,*) ' MATRIX C' 

WRITE (7 , *) ' MATRIX C' 

CALL WRITIT (M,N,C) 

_ 55 continue 

c 

c SET UP MATRIX F = C * C TRANS 

_C 

D0104I=1 , M 
DO104 J=1 , M 
F ( I , J) =0 . 

- D0104K=1,N 

104 F (I , J) =F ( I , J) +C ( I , K) *C(J , K) 
if (ip.lt. 4) goto56 

_ WRITE ( 6 , * ) ' MATRIX F' 

WRITE ( 7 , * ) ' MATRIX F' 

CALL WRITIT (M,M,F) 

56 continue 
~C 

C GET THE INVERSE OF F = FI 

C 

- CALL MATINV (MAX, M, F , MDUM, FI , IOP , DETERM, ISCALE, IPIVOT , IWK) 
WRITE ( 6 , * ) ' DETERM =/ , DETERM , ' ISCALE= / , ISCALE 

WRITE (7 , * ) ' DETERM=' , DETERM , ' ISCALE= / , ISCALE 
_ D0301I=1,M 

DO301J=l , M 
301 FI (I , J) =F ( I , J) 

if (ip. It . 4) goto57 
~ WRITE (6,*)' MATRIX FI' 

WRITE (7 , *) ' MATRIX FI ' 

CALL WRITIT (M,M, FI) 

- 57 continue 

c 

c GET THE PSEUDO INVERSE OF C = CPI = C TRANS * FI 

_C 

DO105IQ=l , N 
DO105J=l,M 
CPI (IQ, J) =0 . 

“ DO105I=l , M 

105 CPI ( IQ , J) =CPI ( IQ , J) +C ( I , IQ) *FI(I, J) 
if ( ip . It . 4 ) goto58 



non non n n non 


WRITE ( 6 , * ) 7 MATRIX CPI 7 
WRITE (7,*)' MATRIX CPI 7 
CALL WRITIT (N, M, CPI ) 

58 continue 

SET UP MATRIX H = PSEUDO INVERSE OF B = BPI 

D0106I=1,M 
D0106J=1 , M 

106 H ( I , J) =BPI ( I , J) 
if (ip. It. 4) goto59 
WRITE(6, *) 7 MATRIX H 7 
WRITE (7 , *) 7 MATRIX H 7 
CALL WRITIT (M,M,H) 

59 continue 

GET INVERSE OF H = HI 

CALL MATINV (MAX ,M,H, MDUM , HI , IOP , DETERM , ISCALE, IPIVOT , IWK) 
WRITE ( 6 , * ) 7 DETERM=' , DETERM, 7 ISCALE= 7 , ISCALE 
WRITE (7 , *) 7 DETERM= 7 ,DETERM, 7 ISCALE= 7 , ISCALE 
D0302I=1 , M 
DO302 J=1 , M 
302 HI (I , J) =H (I , J) 

if ( ip . It . 4 ) goto60 
WRITE (6 , * ) 7 MATRIX HI 7 
WRITE (7,*) 7 MATRIX HI 7 
CALL WRITIT (M,M, HI) 

60 continue 

GET PSEUDO INVERSE OF A = API = CPI * HI * BPI 

D0107IQ=1,N 
DO107 J=1 , M 
API (IQ, J) =0 . 

D0107I=1 , M 
D0107K=1 , M 

107 API (IQ, J) =API ( IQ , J) ++CPI (IQ, I ) *HI ( I , K) *BPI (K, J) 
if ( ip . It . 4 ) goto6 1 

WRITE (6,*) 7 MATRIX API 7 
WRITE (7,*) 7 MATRIX API 7 
CALL WRITIT (N,M, API) 

61 continue 

GET XX = API * Y 

DO108IQ=l,N 
XX (IQ) =0 . 

DO108J=l,M 

108 XX(IQ)=XX(IQ) +API (IQ, J) *Y(J) 

JDUM=1 

if (ip. It . 4) goto62 
WRITE ( 6 , * ) 7 MATRIX XX 7 
WRITE (7 , * ) 7 MATRIX XX 7 
CALL WRITIT (N,JDUM, XX) 

62 continue 

RETURN 

END 

SUBROUTINE WRITIT (MM, NN, XX) 

IMPLICIT REAL* 8 (A-H,0-Z) 



DIMENSION XX (2 1,1) 

1 FORMAT (IX) 

2 FORMAT (10F7. 2) 

WRITE (6,1) 

D0100I=1,MM 

WRITE (6, 2) (XX(I,J) ,J=1,NN) 

WRITE (7 , 2 ) (XX(I,J) ,J=1,NN) 

100 CONTINUE 
RETURN 
END 

SUBROUTINE EXACT ( IP , M, A, Y , B) 

IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION a (250, 136) ,b(136) ,y(250) 

DIMENSION IPIVOT(250) , IWK(250, 2) 

DIMENSION C( 136,1) 

D0100I=1 ,M 

100 C ( I , 1 ) =Y ( I ) 

MAX=250 

MDUM=1 

CALL°MATINV (MAX ,M, A, MDUM , C , IOP , DETERM , ISCALE , IPIVOT , IWK) 
WRITE (6,*) ' DETERM= ' , DETERM, ' ISCALE=' , ISCALE 
WRITE(7,*)' DETERM= ' , DETERM , ' ISCALE=' , ISCALE 
D0101I=1,M 
B (I) =C (1,1) 

101 CONTINUE 

if (ip. It. 4) goto50 

WRITE ( 6 , * ) ' MATRIX B' , (B (I) , 1=1 ,M) 

WRITE ( 7 , * ) ' MATRIX B ' , (B ( I ) , 1=1 , M) 

50 continue 
RETURN 
END 

SUBROUTINE OVER ( IP , M, N , A , Y , B) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION a (250, 136) ,b(136) ,y (250) 

DIMENSION IPIVOT(136) ,IWK(136,2) 

DIMENSION AT A (136,136) ,ATY( 136,1) 

DO200I=l,N 
D0200J=1,N 
ATA(I , J) =0 . 

D0200K=1,M 

200 ATA(I, J)=ATA(I, J)+A(K,I)*A(K, J) 

DO201I=l , N 

ATY (I , 1) =0 . 

D0201K=1,M 

201 ATY (1,1) =ATY (1,1) +A (K, I ) *Y (K) 

MAX=13 6 

MDUM=1 

CALL°MATINV (MAX , N , ATA, MDUM, ATY , IOP , DETERM, ISCALE, IPIVOT , IWK) 
WRITE ( 6 , * ) ' DETERM=' , DETERM, ' ISCALE=' , ISCALE 
WRITE(7,*)' DETERM= ', DETERM , ' I SCALE=' , ISCALE 
D0101I=1,N 
B ( I ) =ATY (1,1) 

101 CONTINUE 

if (ip. It . 4) goto50 

WRITE ( 6 , * ) ' MATRIX B ' , (B ( I) , 1=1 , N) 

WRITE ( 7 , * ) ' MATRIX B ' , (B ( I) , 1=1 , N) 

50 continue 
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RETURN 

END 

subroutine sta t i t ( m , y , yha t , ng , yy , yyha t ) 

— implicit real*8 (a-h,o-z) 
c 

Q *********************************************************** 

_c 

c This subroutine calculates quality of approximation measures 

c this subroutine calculates v, r2 , and vg 

c 

Q *********************************************************** 

c 

dimension y(250) ,yhat(250) 

_ dimension yy (2000) ,yyhat (2000) 

yb=0. 

dol00id=l,m 

yb=yb+y(id) 

100 continue 
yb=yb/ float (m) 
error=0 . 

— dol01id=l,m 

error=error+ (y (id) -yhat(id) ) **2 

101 continue 

_ v=sqrt (error/ float (m) ) /yb*(100. ) 

write (7,*)' error over designs=error = ' , error 
write(7,*)' average y over design = yb =',yb 
write(6,*)' coefficient v (as %)= ' ,v 
write(7,*)' coefficient v (as %)= ',v 
dn=0 . 
dd=0 . 

do7769id=l,m 
dn=dn+(yhat (id) -yb) **2 
dd=dd+(y (id) -yb) **2 
_ 7769 continue 

r2=dn/dd* (100. ) 

write (6,*)' coefficient r2 (as%) = ',r2 
write(7,*)' coefficient r2 (as%) = ',r2 
~ c get vg 

perror=0 . 
yg=0. 

dol55id=l,ng 
yg=yg+yy(id) 

perror=perror+ (yy ( id) -yyhat ( id) ) **2 
155 continue 

yg=yg/float(ng) 

vg=sqrt (perror/f loat (ng) ) /yg* ( 100 . ) 

write(7,*)' sum of residuals squared=perror= ' , perror 
write (7,*)' average y over grid = yg = ',yg 
write (6,*)' coefficient vg = ',vg 
write (7,*)' coefficient vg = ',vg 
return 
end 

SUBROUTINE MATINV (MAX , N, A, M, B, IOP, DETERM, ISCALE, IPIVOT , IWK) MATINV 2 

implicit real*8 (a-h,o-z) 

FI. 3 MATINV 3 

***********************************************************************MATINV 4 

MATINV 5 

PURPOSE - MATINV INVERTS A REAL SQUARE MATRIX A. MATINV 6 

IN ADDITION THE ROUTINE SOLVES THE MATRIX MATINV 7 



c 

c 

c 

-c 

c 

c 

_c 

C USE 

c 

c 

~c 

c 

c 

-C 

c 
c 
_c 
c 
c 
c 
-c 
c 
c 
_c 
c 
c 
c 
~c 
c 
c 
-c 
c 
c 
_c 
c 
c 
c 
_ c 
c 
c 
-c 
c 
c 
_c 
c 
c 
c 
— c 
c 
c 

-C 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 


EQUATION AX=B, WHERE B IS A MATRIX OF CONSTANT 
VECTORS. THERE IS ALSO AN OPTION TO HAVE THE 
DETERMINANT EVALUATED. IF THE INVERSE IS NOT 
NEEDED, USE GELIM TO SOLVE A SYSTEM OF SIMULTANEOUS 
EQUATIONS AND DETFAC TO EVALUATE A DETERMINANT 
FOR SAVING TIME AND STORAGE. 

- CALL MATINV (MAX, N, A,M, B, IOP, DETERM, ISCALE, IPIVOT , IWK) 

MAX - THE MAXIMUM ORDER OF A AS STATED IN THE 

DIMENSION STATEMENT OF THE CALLING PROGRAM. 

N - THE ORDER OF A, 1 . LE.N.LE.MAX. 

A - A TWO-DIMENSIONAL ARRAY OF THE COEFFICIENTS. 
ON RETURN TO THE CALLING PROGRAM, A INVERSE 
IS STORED IN A. 

A MUST BE DIMENSIONED IN THE CALLING PROGRAM 
WITH FIRST DIMENSION MAX AND SECOND DIMENSION 
AT LEAST N. 

M - THE NUMBER OF COLUMN VECTORS IN B. 

M=0 SIGNALS THAT THE SUBROUTINE IS 
USED SOLELY FOR INVERSION, HOWEVER, 

IN THE CALL STATEMENT AN ENTRY CORRE- 
SPONDING TO B MUST BE PRESENT. 

B - A TWO-DIMENSIONAL ARRAY OF THE CONSTANT 
VECTOR B. ON RETURN TO CALLING PROGRAM, 

X IS STORED IN B. B SHOULD HAVE ITS FIRST 
DIMENSION MAX AND ITS SECOND AT LEAST M. 

IOP - COMPUTE DETERMINANT OPTION. 

IOP=0 COMPUTES THE MATRIX INVERSE AND 
DETERMINANT. 

IOP=l COMPUTES THE MATRIX INVERSE ONLY. 

DETERM- FOR IOP=0-IN CONJUNCTION WITH ISCALE 

REPRESENTS THE VALUE OF THE DETERMINANT 
OF A, DET (A) , AS FOLLOWS. 

DET (A) = (DETERM) (10**100 (ISCALE) ) 

THE COMPUTATION DET (A) SHOULD NOT BE 
ATTEMPTED IN THE USER PROGRAM SINCE IF 
THE ORDER OF A IS LARGER AND/OR THE 
MAGNITUDE OF ITS ELEMENTS ARE LARGE (SMALL) , 
THE DET (A) CALCULATION MAY CAUSE OVERFLOW 
(UNDERFLOW) . DETERM SET TO ZERO FOR 
SINGULAR MATRIX CONDITION, FOR EITHER 
I0P=1 , OR 0. SHOULD BE CHECKED BY PROGRAMER 
ON RETURN TO MAIN PROGRAM. 

ISCALE - A SCALE FACTOR COMPUTED BY THE 
SUBROUTINE TO AVOID OVERFLOW OR 
UNDERFLOW IN THE COMPUTATION OF 
THE QUANTITY, DETERM. 

IPIVOT - A ONE DIMENSIONAL INTEGER ARRAY 
USED BY THE SUBPROGRAM TO STORE 
PIVOTOL INFORMATION. IT SHOULD BE 
DIMENSIONED AT LEAST N. IN GENERAL 
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c 


THE USER DOES NOT NEED TO MAKE USE 

MATINV68 

c 


OF THIS ARRAY. 

MATINV69 

c 



MATINV70 

-c 


IWK - A TWO-DIMENSIONAL INTEGER ARRAY OF 

MATINV7 1 

c 


TEMPORARY STORAGE USED BY THE ROUTINE. 

MATINV72 

c 


IWK SHOULD HAVE ITS FIRST DIMENSION 

MATINV73 

c 


MAX, AND ITS SECOND 2. 

MATINV74 

c 



MATINV75 

c 


REQUIRED ROUTINES- 

MATINV76 

c 



MATINV77 

“C 


REFERENCE -FOX,L, AN INTRODUCTION TO NUMERICAL 

MATINV78 

c 


LINEAR ALGEBRA 

MATINV79 

c 



MATINV80 

_c 


STORAGE - 542 OCTAL LOCATIONS 

MATINV81 

c 



MATINV82 

c 


LANGUAGE -FORTRAN 

MATINV83 

c 


LIBRARY FUNCTIONS -ABS 

MATINV84 

c 



MATINV85 

c 


RELEASED - JULY 1973 

MATINV86 

c 



MATINV87 

-c 


LATEST REVISION - JULY 29, 1981 

MATINV88 

c 


COMPUTER SCIECES CORPORATION 

MATINV89 

c 


HAMPTON, VA 

MATINV90 

C***********************************************************************MATINV91 

C 



MATINV92 



DIMENSION IPIVOT(N) ,A(MAX,N) ,B(MAX,N) , IWK (MAX, 2) 

MATINV93 



EQUIVALENCE (IROW,JROW), ( ICOLUM, JCOLUM) , (AMAX, T, SWAP) 

MATINV94 

~C 



MATINV98 

C 


INITIALIZATION 

MATINV99 

C 



MATIN100 




ISCALE=0 

MATIN101 



Rl= ( 10 . 0d+00) **32 

MATIN102 



R2=l. 0d+00/Rl 

MATIN103 



DETERM=1 . 0d+00 

MATIN104 



DO 20 J=1 , N 

MATIN105 



IPIVOT ( J) =0 

MATIN106 

20 

CONTINUE 

MATIN107 

— 


DO 550 1=1, N 

MATIN108 

c 



MATIN109 

c 


SEARCH FOR PIVOT ELEMENT 

MATIN110 

_c 



MATIN111 



AMAX=0 . 0d+00 

MATIN112 



DO 105 J=1 , N 

MATIN113 



IF ( IPIVOT ( J) -1 ) 60, 105, 60 

MATIN114 


60 

DO 100 K=1,N 

MATIN115 



IF (IPIVOT(K) -1) 80, 100, 740 

MATIN116 


80 

TMAX = ABS (A ( J, K) ) 

MATIN117 

— 


IF (AMAX-TMAX) 85,100,100 

MATIN118 


85 

IROW=J 

MATIN119 



ICOLUM=K 

MATIN120 



AMAX=TMAX 

MATIN12 1 


100 

CONTINUE 

MATIN122 


105 

CONTINUE 

MATIN123 



IF (AMAX) 740,106,110 

MATIN124 


106 

DETERM=0 . 0d+00 

MATIN125 



ISCALE=0 

MATIN12 6 



GO TO 740 

MATIN127 

— 

110 

IPIVOT (ICOLUM) = 1 

MATIN128 

c 



MATIN129 

c 


INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

MATIN130 



~C 

140 


200 

210 

_ 250 
260 


C 

C 

-C 

1000 

1010 


1020 


1030 
_ 1040 


1050 

1060 

1070 


1080 


1090 

2000 


2010 

320 
C 

_c 

c 

321 

350 

360 
- 370 

C 

c 


IF ( IROW-ICOLUM) 140, 260, 140 

DETERM=- DETERM 
DO 200 L=1,N 

SWAP=A ( IROW , L) 

A ( IROW , L) =A ( ICOLUM , L) 
A(ICOLUM,L)=SWAP 
CONTINUE 

IF (M) 260, 260, 210 
DO 250 L=1 , M 
SWAP=B ( IROW , L) 

B ( IROW , L) =B ( ICOLUM , L) 

B (ICOLUM, L)=SWAP 
CONTINUE 
IWK (1,1) =IROW 
IWK (1,2) =ICOLUM 
PIVOT=A( ICOLUM, ICOLUM) 

IF(IOP) 740,1000,321 

SCALE THE DETERMINANT 

PIVOTI=PIVOT 

IF (ABS (DETERM) -Rl) 1030,1010,1010 
DETERM=DETERM/ Rl 
ISCALE=ISCALE+1 

IF (ABS (DETERM) -Rl) 1060,1020,1020 
DETERM=DETERM/R1 
ISCALE=ISCALE+1 
GO TO 1060 

IF (ABS (DETERM) -R2) 1040,1040,1060 

DETERM=DETERM*R1 

ISCALE=I SCALE- 1 

IF (ABS (DETERM) -R2) 1050, 1050, 1060 

DETERM=DETERM*R1 

ISCALE=I SCALE- 1 

IF(ABS (PIVOTI) -Rl) 1090,1070,1070 

PIVOTI=PIVOTI/Rl 

ISCALE=ISCALE+1 

IF (ABS (PIVOTI) -Rl) 320,1080,1080 
PIVOTI=PIVOTI /Rl 
ISCALE=ISCALE+1 
GO TO 320 

IF (ABS (PIVOTI) -R2) 2000, 2000, 320 
PIVOTI=PIVOTI *R1 
ISCALE=ISCALE-1 

IF (ABS (PIVOTI ) -R2) 2010, 2010, 320 
PIVOTI=PI VOTI *R1 
ISCALE=ISCALE-1 
DETERM=DETERM*PIVOTI 

DIVIDE PIVOT ROW BY PIVOT ELEMENT 

A (ICOLUM, ICOLUM) =1. 0d+00 

DO 350 L=1,N 

A ( ICOLUM, L)=A( ICOLUM, L) /PIVOT 
IF (M) 380, 380, 360 
DO 370 L=1,M 

B ( ICOLUM , L) =B ( ICOLUM , L) /PIVOT 
REDUCE NON-PIVOT ROWS 
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c 


c 

c 

_C 


~c 

c 

c 

-C 

c 
c 
_c 
c 
c 
c 
— c 
c 
c 

-C 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 

_c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 


380 

DO 550 Ll=l , N 

IF (Ll-ICOLUM) 400, 550, 400 

400 

T=A (LI , ICOLUM) 

A (LI , ICOLUM) =0 . 0d+00 
DO 450 L=1 , N 


450 

A (LI , L) =A (LI , L) -A (ICOLUM, 
IF (M) 550, 550, 460 

L) *T 

460 

DO 500 L=1 , M 


500 

B(L1,L)=B(L1,L) -B( ICOLUM, 

L) *T 

550 

CONTINUE 



INTERCHANGE COLUMNS 

DO 710 1=1, N 
L=N+1-I 

IF (IWK(L, 1)-IWK(L, 2) ) 630,710,630 
630 JROW=IWK ( L , 1 ) 

JCOLUM=IWK ( L , 2 ) 

DO 705 K=1,N 
SWAP=A(K, JROW) 

A (K , JROW) =A (K, JCOLUM) 

A ( K , J COLUM ) =SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 

ROUTINE NAME - HC318=EPSLON 
FROM EISPACK 


LATEST REVISION 

- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 

PURPOSE 

- THE FORTRAN FUNCTION EPSLON ESTIMATES UNIT 
ROUNDOFF IN QUANTITIES OF SIZE X. 

USAGE 

- VARIABLE = EPSLON (X) 

ARGUMENTS X 

- IS A REAL INPUT VARIABLE WHICH REPRESENTS 
QUANTITIES OF SIZE IN WHICH UNIT ROUNDOFF 
WILL BE ESTIMATED. 

REQUIRED ROUTINES 

- NONE 


REMARKS 1. IT SHOULD BE NOTED THAT EPSLON IS A FUNCTION 
DESIGNED TO BE CALLED BY ROUTINES IN THE 
EISPACK VERSION 3. 


THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL 
SYSTEMS SATISFYING THE FOLLOWING TWO 
ASSUMPTIONS, 

A. THE BASE USED IN REPRESENTING FLOATING 
POINT NUMBERS IS NOT A POWER OF THREE. 

B. THE QUANTITY A IN STATEMENT 10 IS 
REPRESENTED TO THE ACCURACY USED IN FLOATING 
POINT VARIABLES THAT ARE STORED IN MEMORY. 


MATIN191 
MATIN192 
MATIN193 
MATIN194 
MATIN195 
MATIN196 
MATIN197 
MATIN198 
MATIN199 
MATIN200 
MATIN201 
MATIN202 
MATIN203 
MATIN204 
MATIN205 
MATIN206 
MATIN207 
MATIN208 
MATIN209 
MATIN210 
MATIN2 11 
MATIN2 12 
MATIN2 13 
MATIN214 
MATIN215 
MATIN216 
MATIN2 17 
EPSLON 2 
EPSLON 3 
EPSLON 4 
EPSLON 5 
EPSLON 6 
EPSLON 7 
EPSLON 8 
EPSLON 9 
EPSLONIO 
EPSLON11 
EPSLON12 
EPSLON 13 
EPSLON 14 
THE EPSLON 15 
EPSLON 16 
EPSLON 17 
EPSLON 18 
EPSLON19 
EPSLON20 
EPSLON21 
EPSLON22 
EPSLON23 
EPSLON24 
EPSLON25 
EPSLON26 
EPSLON27 
EPSLON28 
EPSLON29 
EPSLON30 
EPSLON31 
EPSLON32 
EPSLON33 
EPSLON34 



c 

C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE 

C INTENDED TO FORCE OPTIMIZING COMPILERS TO 

C GENERATE CODE SATISFYING ASSUMPTION 2. 

C 

C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE 

C THAT, 

C 

C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, 

C 

C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, 

C 

C C IS NOT EXACTLY EQUAL TO ONE, 

C 

C EPS MEASURES THE SEPARATION OF 1.0 FROM THE 

C NEXT LARGER FLOATING POINT NUMBER. 

C 

C EXAMPLE : 

C PROGRAM TR( OUTPUT, TAPE6=OUTPUT) 

C REAL X 

C X = 4. 

C A = EPSLON (X) 

C WRITE (6,100) A 

C100 FORMAT (5H0A = ,G22.14) 

C STOP 

C END 

C OUTPUT : 

CA = . 56843418860808E-13 

C 

C 

C*F45V1P0* 

REAL* 8 FUNCTION EPSLON (X) 

C 

REAL* 8 X 

REAL* 8 A, B, C, EPS 
A = 4.0E0/3.0E0 
10 B = A - 1.0E0 
C = B + B + B 
EPS = ABS(C-l.OEO) 

IF (EPS .EQ. 0.0E0) GO TO 10 
EPSLON = EPS*ABS (X) 

RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


ROUTINE NAME 
FROM El SPACE 


LATEST REVISION 


PURPOSE 


PF260=QZHES 


AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 


THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
GENERAL MATRICES AND REDUCES ONE OF THEM TO 
UPPER HESSENBERG FORM AND THE OTHER TO UPPER 
TRIANGULAR FORM USING ORTHOGONAL 
TRANSFORMATIONS. IT IS USUALLY FOLLOWED BY 
QZIT (PF261) , QZVAL(PF262) AND, POSSIBLY, 
QZVEC (PF263 ) . 


EPSLON35 
EPSLON36 
EPSLON37 
EPSLON38 
EPSLON39 
EPSLON40 
EPSLON41 
EPSLON42 
EPSLON43 
EPSLON44 
EPSLON45 
EPSLON46 
EPSLON47 
EPSLON48 
EPSLON49 
EPSLON50 
EPSLON51 
EPSLON52 
EPSLON53 
EPSLON54 
EPSLON55 
EPSLON56 
EPSLON57 
EPSLON58 
EPSLON59 
EPSLON60 
EPSLON61 
EPSLON 6 2 
EPSLON 6 3 
EPSLON64 
EPSLON 6 5 
EISPAK 
EISPAK32 
EISPAK 
EISPAK 
EISPAK35 
EISPAK36 
EISPAK37 
EISPAK38 
EISPAK39 
EISPAK40 
EISPAK41 
EISPAK42 
EISPAK43 
QZHES 2 
QZHES 3 
QZHES 4 
QZHES 5 
QZHES 6 
QZHES 7 
QZHES 8 
QZHES 9 
QZHES 10 
QZHES 11 
QZHES 12 
QZHES 13 
QZHES 14 
QZHES 15 
QZHES 16 
QZHES 17 



c 

c 

c 

c 

c 

c 

_c 

c 


"C 

c 


c 

c 

n 

w 

-c 

c 


c 

c 


~c 

c 

c 


c 

-C 

c 


c 

r* 

"C 

c 


USAGE 


ARGUMENTS 


NM 


N 


- CALL QZHES ( NM, N, A, B, MATZ, Z) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT . 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL GENERAL MATRIX. 
MUST BE OF DIMENSION NM X N. 


QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 

QZHES 


B 


MATZ 


ON OUTPUT A HAS BEEN REDUCED TO UPPER 
HESSENBERG FORM. THE ELEMENTS BELOW THE FIRSTQZHES 
SUBDIAGONAL HAVE BEEN SET TO ZERO. QZHES 

QZHES 

ON INPUT B CONTAINS A REAL GENERAL MATRIX. QZHES 
MUST BE OF DIMENSION NM X N. QZHES 

QZHES 

ON OUTPUT B HAS BEEN REDUCED TO UPPER QZHES 

TRIANGULAR FORM. THE ELEMENTS BELOW THE MAIN QZHES 40 
DIAGONAL HAVE BEEN SET TO ZERO. QZHES 41 

QZHES 42 

ON INPUT MATZ SHOULD BE SET TO .TRUE. IF THE QZHES 43 


18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 


RIGHT HAND TRANSFORMATIONS ARE TO BE 
ACCUMULATED FOR LATER USE IN COMPUTING 
EIGENVECTORS, AND TO .FALSE. OTHERWISE. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE RIGHT 
HAND TRANSFORMATIONS IF MATZ HAS BEEN SET TO 
.TRUE. OTHERWISE, Z IS NOT REFERENCED. 

MUST BE OF DIMENSION NM X N. 


REQUIRED ROUTINES - NONE 


REMARKS 1. 


THIS SUBROUTINE IS THE FIRST STEP OF THE QZ 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL. 10, 
241-256(1973) BY MOLER AND STEWART. 


EXAMPLE : 

PROGRAM TQZHES (OUTPUT, TAPE6=OUTPUT) 
DIMENSION A ( 5 , 5 ) , Z ( 5 , 5 ) ,B(5,5) 
LOGICAL MATZ 

N = 5 
NM = 5 

MATZ = .TRUE. 


DATA A /10. ,2. ,3. ,2*1. ,2 
* 1. ,-l. ,1. ,2. ,1. ,9 


, 12 . , 1 . , 2 • , 1 . ^ 3 • , 1 * , 1 1 • , 
,3*1. ,-l. ,1. ,15. / 


DATA B /12. ,1. ,-l. ,2. ,2*1. ,14. ,1. ,-l. ,1. ,“1. ,1. , 

* 16. ,-l. , 1. , 2. ,-l. ,-l. , 12. ,-l. ,3*1. ,-l. ,11. 


/ 


CALL QZHES (NM, N, A, B , MATZ , Z) 

WRITE (6, 100) ( (A ( I , J) ,1=1,5) , J=l,5) , ( (B(I,J) ,1=1,5) ,J=1,5) , 


QZHES 44 
QZHES 45 
QZHES 46 
QZHES 47 
QZHES 48 
QZHES 49 
QZHES 50 
QZHES 51 
QZHES 52 
QZHES 53 
QZHES 54 
QZHES 55 
QZHES 56 
QZHES 57 
QZHES 58 
QZHES 59 
QZHES 60 
QZHES 61 
QZHES 62 
QZHES 63 
QZHES 64 
QZHES 
QZHES 
QZHES 
QZHES 
QZHES 69 
QZHES 70 
QZHES 71 
QZHES 72 
QZHES 73 
QZHES 
QZHES 
QZHES 
QZHES 


65 

66 

67 

68 


74 

75 

76 

77 



c * ( (Z(I,J) ,1=1,5) ,J=1,5) 

C100 FORMAT (1H , 5H A = /5(1H , 5 (G8 . 2 , 2X) / ) ) 
C * 5H B = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) 

- C * 5H Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 


c 

STOP 





c 

END 





_c 






c 

OUTPUT : 





c 






c 

A = 





~c 

-9.9 

4.1 

0. 

0. 

0. 

c 

-2.4 

11. 

-3.0 

0. 

0. 


.91 

.26 

-13 . 

3.3 

0. 

-C 

-3.8 

2 . 0 

1.7 

-11. 

2 . 6 

c 

2.7 

-1.5 

-.99 

1.4 

-11. 

c 

B = 






-12. 

0. 

0. 

0. 

0. 

~c 

2.3 

16. 

0. 

0. 

0. 

c 

-.34 

-3.0 

-12. 

0. 

0. 

c 

-3.8 

.80 

-1.5 

-10. 

0 . 

-c 

2.5 

-1.4 

-1.5 

-1.5 

-13. 

c 

Z = 





r* 

1.0 

0. 

0. 

0. 

0. 

<-» 

0. 

.26 

.95 

-.14 

- . 70E-01 



c 

0. 

. 87E-01 

- . 24E-01 

.43 

-.90 

c 

0. 

. 24E-01 

. 16 

.89 

.43 

c 

0. 

-.96 

.26 

. 22E-01 

- . 89E-01 


C 



SUBROUTINE QZHES (NM, N , A, B , MATZ , Z) 


QZHES 78 
QZHES 79 
QZHES 80 
QZHES 81 
QZHES 82 
QZHES 83 
QZHES 84 
QZHES 85 
QZHES 86 
QZHES 87 
QZHES 88 
QZHES 89 
QZHES 90 
QZHES 91 
QZHES 92 
QZHES 93 
QZHES 94 
QZHES 95 
QZHES 96 
QZHES 97 
QZHES 98 
QZHES 99 
QZHES100 
QZHES101 
QZHES102 
QZHES103 
QZHES104 
QZHES105 
QZHES106 
EISP6685 


implicit real*8 (a-h,o-z) 

INTEGER I , J , K , L , N , LB , LI , NM , NK1 , NM1 , NM2 
REAL*8 A(NM,N) , B (NM, N) , Z (NM,N) 

REAL* 8 R, S , T , U1 , U2 , VI , V2 , RHO 
LOGICAL MATZ 

IF (.NOT. MATZ) GO TO 10 
“C 

DO 3 J = 1, N 
C 

- DO 2 I = 1, N 

Z(I,J) = 0.0E0 

2 CONTINUE 

Z(J,J) = 1.0E0 

3 CONTINUE 

q REDUCE B TO UPPER TRIANGULAR FORM 

- 10 IF (N .LE. 1) GO TO 170 

NM1 = N - 1 

i 

_ DO 100 L = 1, NM1 

LI = L + 1 
S = 0.0E0 

n 

DO 20 I = LI, N 

S = S + ABS (B ( I , L) ) 

20 CONTINUE 
~C 

IF (S .EQ. 0.0E0) GO TO 100 
S = S + ABS (B (L, L) ) 


EISP6686 
EISP6687 
EISP66 
EISP66 
EISP6690 
EISP6691 
EISP6692 
EISP6693 
EISP6694 
EISP6695 
EISP6696 
EISP6697 
EISP6698 
EISP6699 
EISP6700 
EISP6701 
EISP6702 
EISP6703 
EISP67 04 
EISP67 05 
EISP67 06 
EISP6707 
EISP6708 
EISP6709 
EISP67 10 
EISP67 11 
EISP6712 
EISP67 13 
EISP67 14 


o o 


R = O.OEO EISP67 15 

C EISP67 16 

DO 25 I = L, N EISP6717 

B ( I , L) = B ( I , L) / S EISP67 18 

R = R + B ( I , L) **2 EISP67 19 

25 CONTINUE EISP6720 

C EISP6721 

R = SIGN (SQRT (R) , B (L, L) ) EISP6722 

B (L, L) = B (L, L) + R EISP6723 

RHO = R * B (L, L) EISP6724 

C EISP6725 

DO 50 J = LI, N EISP6726 

T = O.OEO EISP6727 

C EISP6728 

DO 30 I = L, N EISP6729 

T = T + B (I , L) * B ( I , J) EISP6730 

30 CONTINUE EISP6731 

C EISP6732 

T = -T / RHO EISP673 3 

C EISP6734 

DO 40 I = L, N EISP6735 

B ( I , J) = B ( I , J) + T * B ( I , L) EISP6736 

40 CONTINUE EISP6737 

C EISP6738 

50 CONTINUE EISP6739 

C EISP6740 

DO 80 J = 1, N EISP6741 

T = O.OEO EISP6742 

C EISP6743 

DO 60 I = L, N EISP6744 

T = T + B ( I , L) * A ( I , J) EISP6745 

60 CONTINUE EISP6746 

C EISP6747 

T = -T / RHO EISP6748 

C EISP6749 

DO 70 I = L, N EISP6750 

A ( I , J) = A ( I , J) + T * B (I , L) EISP6751 

70 CONTINUE EISP6752 

C EISP6753 

80 CONTINUE EISP6754 

C EISP6755 

B (L, L) = -S * R EISP6756 

C EISP6757 

DO 90 I = LI, N EISP6758 

B (I , L) = O.OEO EISP6759 

90 CONTINUE EISP6760 

C EISP6761 

100 CONTINUE EISP6762 

REDUCE A TO UPPER HESSENBERG FORM, WHILE EISP6763 

KEEPING B TRIANGULAR EISP6764 

IF (N .EQ. 2) GO TO 170 EISP6765 

NM2 = N - 2 EISP6766 

C EISP67 67 

DO 160 K = 1, NM2 EISP6768 

NK1 = NM1 - K EISP6769 

C FOR L=N-1 STEP -1 UNTIL K+l DO — EISP6770 

DO 150 LB = 1, NK1 EISP6771 

L = N - LB EISP6772 

LI = L + 1 EISP6773 

C ZERO A (L+l , K) EISP6774 





S = ABS (A (L, K) ) + ABS (A (LI , K) ) 

EISP6775 



IF (S .EQ. 0.0E0) GO TO 150 

EISP6776 



U1 = A(L,K) / S 

EISP6777 

— 


U2 = A (LI , K) / S 

EISP6778 



R = SIGN ( SQRT (U1*U1+U2*U2 ) ,U1) 

EISP6779 



VI = -(Ul + R) / R 

EISP6780 



V2 = -U2 / R 

EISP6781 



U2 = V2 / VI 

EISP6782 

c 



EISP6783 



DO 110 J = K, N 

EISP6784 



T = A (L, J) + U2 * A (LI , J) 

EISP6785 



A (L, J) = A (L, J) + T * VI 

EISP6786 



A (LI , J) = A (LI , J) + T * V2 

EISP6787 

■ — 

110 

CONTINUE 

EISP6788 

c 



EISP6789 



A (LI , K) = 0.0E0 

EISP6790 

c 



EISP6791 



DO 120 J = L, N 

EISP6792 



T = B (L, J) + U2 * B (LI , J) 

EISP6793 



B (L, J) = B (L, J) + T * VI 

EISP6794 

— 


B (LI , J) = B (LI , J) + T * V2 

EISP6795 


120 

CONTINUE 

EISP6796 

c 


ZERO B (L+l , L) 

EISP6797 




S = ABS (B (LI , LI) ) + ABS (B (LI , L) ) 

EISP6798 



IF (S .EQ. 0.0E0) GO TO 150 

EISP6799 



Ul = B (LI , LI) / S 

EISP6800 



U2 = B (LI , L) / S 

EISP6801 



R = SIGN ( SQRT (U1*U1+U2 *U2 ) ,U1) 

EISP6802 



VI = -(Ul + R) / R 

EISP6803 



V2 = -U2 / R 

EISP6804 

— 


U2 = V2 / VI 

EISP6805 

c 



EISP6806 



DO 130 1=1, LI 

EISP6807 



T = B ( I , LI ) + U2 * B (I , L) 

EISP6808 



B (I , LI) = B (I , LI) + T * VI 

EISP6809 



B (I , L) = B ( I , L) + T * V2 

EISP6810 


130 

CONTINUE 

EISP6811 

-C 



EISP6812 



B (LI , L) = 0.0E0 

EISP6813 

c 



EISP6814 

— 


DO 140 I = 1, N 

EISP6815 



T = A(I , LI) + U2 * A ( I , L) 

EISP6816 



A ( I , LI) = A(I , LI) + T * VI 

EISP6817 



A ( I , L) = A (I , L) + T * V2 

EISP6818 


140 

CONTINUE 

EISP6819 

c 



EISP6820 



IF (.NOT. MATZ) GO TO 150 

EISP682 1 

~c 



EISP6822 



DO 145 I = 1, N 

EISP6823 



T = Z ( I , LI ) + U2 * Z (I , L) 

EISP6824 

« 


Z (I , LI) = Z (I , LI) + T * VI 

EISP6825 



Z ( I , L) = Z ( I , L) + T * V2 

EISP6826 


145 

CONTINUE 

EISP6827 

_c 



EISP6828 


150 

CONTINUE 

EISP6829 

c 



EISP6830 


160 

CONTINUE 

EISP683 1 

-c 



EISP6832 


170 

RETURN 

EISP6833 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

EISP6834 



c 
c 
-c 
c- 
c 
_c 
c 
c 
c 
~ c 
c 
c 

-C 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 



END 

ROUTINE NAME 
FROM EISPACK 


LATEST REVISION 


PURPOSE 


USAGE 

ARGUMENTS NM 


N 

A 


B 


EPS1 


- PF261=QZIT 


- AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 

- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN UPPER HESSENBERG 
FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 
IT REDUCES THE HESSENBERG MATRIX TO 

QUASI -TRIANGULAR FORM USING ORTHOGONAL 
TRANSFORMATIONS WHILE MAINTAINING THE 
TRIANGULAR FORM OF THE OTHER MATRIX. IT IS 
USUALLY PRECEDED QZHES(PF260) AND FOLLOWED 
BY QZVAL (PF262 ) AND, POSSIBLY, QZVEC (PF263 ) . 

- CALL QZIT (NM , N , A, B , EPS1 , MATZ , Z , IERR) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT . 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL UPPER HESSENBERG 
MATRIX. 

MUST BE OF DIMENSION NM X N. 


EISP6835 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 


ON OUTPUT A HAS BEEN REDUCED TO 
OUASI -TRIANGULAR FORM. THE ELEMENTS BELOW THEQZIT 

_ a \TOi mrjA /’\7T r P 


FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO 
CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. 

ON INPUT B CONTAINS A REAL UPPER TRIANGULAR 
MATRIX. 

MUST BE OF DIMENSION NM X N. 


QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 

ON OUTPUT B IS STILL IN UPPER TRIANGULAR QZIT 
FORM, ALTHOUGH ITS ELEMENTS HAVE BEEN ALTERED. QZIT 
THE LOCATION B(N,1) IS USED TO STORE EPS1 QZIT 
TIMES THE NORM OF B FOR LATER USE BY QZVAL QZIT 
QZVAL (PF2 62) AND QZVEC (PF263) . QZIT 

ON INPUT EPS1 IS A TOLERANCE USED TO DETERMINEQZIT 
NEGLIGIBLE ELEMENTS. EPS1 = 0.0 (OR NEGATIVE) QZIT 
MAY BE INPUT, IN WHICH CASE AN ELEMENT WILL BEQZIT 
NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF QZIT 
ERROR TIMES THE NORM OF ITS MATRIX. IF THE QZIT 
INPUT EPS1 IS POSITIVE, THEN AN ELEMENT WILL QZIT 
BE CONSIDERED NEGLIGIBLE IF IT IS LESS THAN QZIT 
EPS1 TIMES THE NORM OF ITS MATRIX. A POSITIVEQZIT 
VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, QZIT 
BUT LESS ACCURATE RESULTS. Q ZIT 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


u o 


c 

c 

c 

-c 

c 

c 

_c 

c 

c 

c 

~c 

c 

c 

-c 

c 

€ 

_c 

c 

c 

c 

-c 

c 

c 

_c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 

_c 

c 

c 

c 

~c 

c 

c 

-c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 

-C 

C99 


C100 

c 


MATZ 


QZIT 

ON INPUT MATZ SHOULD BE SET TO .TRUE. IF THE QZIT 
RIGHT HAND TRANSFORMATIONS ARE TO BE QZIT 
ACCUMULATED FOR LATER USE IN COMPUTING QZIT 
EIGENVECTORS, AND TO .FALSE. OTHERWISE. QZIT 

QZIT 

ON INPUT Z CONTAINS, IF MATZ HAS BEEN SET TO QZIT 
•TRUE., THE TRANSFORMATION MATRIX PRODUCED IN QZIT 
THE REDUCTION BY QZHES (PF260) , IF PERFORMED, QZIT 
OR ELSE THE IDENTITY MATRIX. IF MATZ HAS BEENQZIT 


SET TO .FALSE., Z IS NOT REFERENCED. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE 
RIGHT HAND TRANSFORMATIONS (FOR BOTH STEPS) 
MATZ HAS BEEN SET TO .TRUE.. 


I ERR - 


QZIT 

QZIT 

QZIT 

QZIT 

IFQZIT 

QZIT 

QZIT 

QZIT 

QZIT 


ON OUTPUT I ERR IS SET TO 
ZERO FOR NORMAL RETURN. 

J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED QZIT 
WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 


REQUIRED ROUTINES - HC318=EPSLON 


REMARKS 


THIS SUBROUTINE IS THE SECOND STEP OF THE QZ 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL. 10, 
241-256(1973) BY MOLER AND STEWART, AS 
MODIFIED IN TECHNICAL NOTE NASA TN 
D-7 3 05 ( 1973 ) BY WARD. 


EXAMPLE : 

PROGRAM TQZIT (OUTPUT, TAPE6=OUTPUT) 
DIMENSION A (5 , 5) ,B(5,5) ,Z(5,5) 
LOGICAL MATZ 


N = 5 
NM = 5 
MATZ = 
EPS1 = 


.TRUE. 

0.0E0 


DATA A /10. ,2. ,3. ,2*1. 
k 1. , -1. , 1. , 2 . , 1. 


,2. ,12. ,1. ,2. ,1. ,3. ,1. ,11. 
,9. ,3*1. ,-l. ,1. ,15. / 


C 

C 

C 

c 


DATA B /12. ,1. ,-l. ,2. ,2*1. ,14. ,1. ,-l. ,1. ,-l. ,1. , 

* 16. , — 1 . ,1. ,2. ,-l. ,~1* ,12. ,-l. ,3*1. , -1 . ,11* / 

CALL QZHES (NM, N , A, B , MATZ , Z) 

CALL QZIT (NM, N, A, B , EPS1 , MATZ , Z , I ERR) 

WRITE (6, 99) I ERR 
FORMAT (1H1,8H I ERR = ,14) 

WRITE (6, 100) ( (A (I, J) ,1=1,5) ,J=1,5) , ( (B(I,J) ,1=1,5) , J=l,5) , 

* ((Z(I,J) ,1=1,5) ,J=1,5) 

FORMAT (1H ,5H A = /5(1H , 5 (G8 . 2 , 2X) / ) ) 

* 5H B = /5(1H , 5 (G8 . 2 , 2X) / ) 

* 5H Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 

STOP 

END 


61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 


QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 
QZIT 100 
QZIT 101 
QZIT 102 
QZIT 103 
QZIT 104 
QZIT 105 
QZIT 106 
QZIT 107 
QZIT 108 
QZIT 109 
QZIT 110 
QZIT 111 
QZIT 112 
QZIT 113 
QZIT 114 
QZIT 115 
QZIT 116 
QZIT 117 
QZIT 118 
QZIT 119 
QZIT 120 



OUTPUT : 


C 

C 


QZIT 121 
QZIT 122 


C 

I ERR = 

0 




QZIT 

123 

-c 

A = 





QZIT 

124 

c 

-15. 

-1.3 

0 . 

0 . 

0 . 

QZIT 

125 

c 

1.1 

7.4 

0 . 

0 . 

0 . 

QZIT 

126 

c 

1.5 

-1.5 

-16. 

0 . 

0 . 

QZIT 

127 

c 

-2.2 

.96 

1.0 

-10. 

0. 

QZIT 

128 


-2 . 6 

-.31 

1.2 

1.7 

-8.6 

QZIT 

129 

c 

B = 





QZIT 

130 

“C 

-9.9 

0 . 

0 . 

0 . 

. 31E-12 

QZIT 

131 

c 

-.29 

17. 

0 . 

0 . 

0. 

QZIT 

132 

c 

1.3 

-2.1 

-14. 

0 . 

0. 

QZIT 

133 

c 

-1.9 

1.7 

.96 

-11. 

0. 

QZIT 

134 

c 

-2.6 

-.32 

1.3 

2.1 

-13. 

QZIT 

135 

'C 

Z = 





QZIT 

136 

c 

.28 

- . 7 IE-01 

.16 

-.24 

-.91 

QZIT 

137 

c 

.52 

-.24 

-.66 

.48 

- . 64E-01 

QZIT 

138 

c 

.49 

.56 

.49 

.45 

. 75E-01 

QZIT 

139 

c 

-.60 

.48 

-.29 

.44 

-.38 

QZIT 

140 

~c 

-.25 

-.63 

.45 

.57 

- . 94E-01 

QZIT 

141 

c 






QZIT 

142 








143 

c — 

SUBROUTINE QZIT (NM 

, N , A, B , 

EPS 1, MATZ, 

Z , IERR) 

EISP6836 


C 

implicit real*8 (a-h,o-z) 

INTEGER I , J, K, L, N, EN, K1 , K2 , LD, LL, LI , NA,NM, ISH, ITN, ITS, KM1 , LM1 , 
X ENM2 , I ERR, LORI , ENORN 

REAL* 8 A(NM,N) ,B(NM,N) ,Z(NM,N) 

REAL* 8 R, S , T, A1 , A2 , A3 , EP, SH,U1 , U2 , U3 , VI, V2 , V3 , ANI , All , 

X A12 , A2 1 , A22 , A3 3 , A3 4 , A43 , A44 , BNI , Bll , B12 , B22 , B3 3 , B34 , 

X B4 4 , EPSA , EPSB , EPS1 , ANORM , BNORM , EPSLON 

LOGICAL MATZ , NOTLAS 
I ERR = 0 

C COMPUTE EPSA, EPSB 

ANORM = 0.0E0 
BNORM = 0.0E0 

"C 

DO 30 I = 1, N 
ANI = 0.0E0 

- IF (I .NE. 1) ANI = ABS (A(I,I-1) ) 

BNI = 0.0E0 

C 

_ DO 20 J = I, N 

ANI = ANI + ABS (A ( I , J) ) 

BNI = BNI + ABS (B ( I , J) ) 

20 CONTINUE 
~C 

IF (ANI .GT. ANORM) ANORM = ANI 
IF (BNI .GT. BNORM) BNORM = BNI 

- 30 CONTINUE 
C 

IF (ANORM .EQ. 0.0E0) ANORM = 1.0E0 
___ IF (BNORM .EQ. 0.0E0) BNORM = 1.0E0 

EP = EPS1 

IF (EP .GT. 0.0E0) GO TO 50 


C USE ROUNDOFF LEVEL IF EPS1 IS ZERO 

EP = EPSLON (1.0E0) 

50 EPSA = EP * ANORM 
EPSB = EP * BNORM 


EISP6837 
EISP6838 
EISP6839 
EISP68 
EISP68 
EISP6842 
EISP6843 
EISP6844 
EISP6845 
EISP6846 
EISP6847 
EISP6848 
EISP6849 
EISP6850 
EISP6851 
EISP6852 
EISP6853 
EISP6854 
EISP6855 
EISP6856 
EISP6857 
EISP6858 
EISP6859 
EISP6860 
EISP6861 
EISP6862 
EISP6863 
EISP6864 
EISP6865 
EISP6866 
EISP6867 
EISP6868 
EISP6869 
EISP6870 
EISP687 1 



o o 




REDUCE A TO QUASI -TRIANGULAR FORM, WHILE EISP6872 

KEEPING B TRIANGULAR EISP6873 

LORI = 1 EISP6874 

ENORN = N EISP6875 

EN = N EISP6876 

ITN = 30*N EISP6877 

BEGIN QZ STEP EISP6878 

60 IF (EN .LE. 2) GO TO 1001 EISP6879 

IF (.NOT. MATZ) ENORN = EN EISP6880 

ITS = 0 EISP6881 

NA = EN - 1 EISP6882 

ENM2 = NA - 1 EISP6883 

70 ISH = 2 EISP6884 

CHECK FOR CONVERGENCE OR REDUCIBILITY . EISP6885 

FOR L=EN STEP -1 UNTIL 1 DO — EISP6886 

DO 80 LL = 1, EN EISP6887 

LM1 = EN - LL EISP6888 

L = LM1 + 1 EISP6889 

IF (L .EQ. 1) GO TO 95 EISP6890 

IF (ABS (A (L, LM1) ) .LE. EPSA) GO TO 90 EISP6891 

80 CONTINUE EISP6892 

EISP6893 

90 A(L, LM1) = 0.0E0 EISP6894 

IF (L .LT. NA) GO TO 95 EISP6895 

1-BY-l OR 2-BY-2 BLOCK ISOLATED EISP6896 

EN = LM1 EISP6897 

GO TO 60 EISP6898 

CHECK FOR SMALL TOP OF B EISP6899 

95 LD = L EISP6900 

100 LI = L + 1 EISP6901 

Bll = B(L,L) EISP6902 

IF (ABS (Bll) .GT. EPSB) GO TO 120 EISP6903 

B(L,L) = 0.0E0 EISP6904 

S = ABS (A(L, L) ) + ABS (A (LI , L) ) EISP6905 

U1 = A(L,L) / S EISP6906 

U2 = A (LI , L) / S EISP6907 

R = SIGN (SQRT (U1*U1+U2*U2 ) ,U1) EISP6908 

VI = -(Ul + R) / R EISP6909 

V2 = -U2 / R EISP6910 

U2 = V2 / VI EISP6911 

EISP6912 

DO 110 J = L, ENORN EISP6913 

T = A (L, J) + U2 * A (LI , J) EISP6914 

A (L, J) = A (L, J) + T * VI EISP6915 

A (LI , J) = A (LI , J) + T * V2 EISP6916 

T = B (L, J) + U2 * B (LI , J) EISP6917 

B (L, J) = B (L, J) + T * VI EISP6918 

B (LI , J) = B (LI , J) + T * V2 EISP6919 

110 CONTINUE EISP6920 

EISP6921 

IF (L .NE. 1) A (L, LM1) = -A (L, LM1) EISP6922 

LM1 = L EISP6923 

L = LI EISP6924 

GO TO 90 EISP6925 

120 All = A (L, L) / Bll EISP6926 

A21 = A (LI , L) / Bll EISP6927 

IF (ISH .EQ. 1) GO TO 140 EISP6928 

ITERATION STRATEGY EISP6929 

IF (ITN .EQ. 0) GO TO 1000 EISP6930 

IF (ITS .EQ. 10) GO TO 155 EISP6931 



non 


C DETERMINE TYPE OF SHIFT 

B22 = B (LI , LI) 

IF (ABS (B22 ) .LT. EPSB) B22 = EPSB 
B33 = B (NA, NA) 

IF (ABS (B33 ) .LT. EPSB) B33 = EPSB 
B44 = B ( EN , EN ) 

IF (ABS (B44 ) .LT. EPSB) B44 = EPSB 

A3 3 = A (NA, NA) / B33 

A3 4 = A (NA, EN) / B44 

A43 = A(EN,NA) / B33 

A44 = A ( EN , EN ) / B44 

B34 = B (NA, EN) / B44 

T = 0.5E0 * (A43 * B34 - A33 - A44) 
R = T * T + A34 * A43 - A33 * A44 
IF (R -LT. 0.0E0) GO TO 150 


C DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A 

ISH = 1 
R = SQRT(R) 

SH = -T + R 
S = -T - R 


IF (ABS ( S-A44 ) .LT. ABS(SH-A44)) SH = S 

LOOK FOR TWO CONSECUTIVE SMALL 

SUB-DIAGONAL ELEMENTS OF A. 

FOR L=EN-2 STEP -1 UNTIL LD DO — 

DO 130 LL = LD, ENM2 
L = ENM2 + LD - LL 
IF (L .EQ. LD) GO TO 140 
LM1 = L - 1 
LI = L + 1 
T = A (L, L) 

IF (ABS (B (L, L) ) .GT. EPSB) T = T - SH * B(L,L) 

IF (ABS (A (L, LM1) ) .LE. ABS (T/A (LI , L) ) * EPSA) GO TO 100 
130 CONTINUE 

140 A1 = All - SH 
A2 = A21 

IF (L .NE. LD) A (L, LM1) = -A (L, LM1) 

GO TO 160 

DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A 

150 A12 = A(L, LI) / B22 
A22 = A (LI , LI) / B22 
B12 = B (L, LI) / B22 

A1 = ( (A3 3 - All) * (A44 - All) - A34 * A43 + A43 * B34 * All) 

X / A21 + A12 - All * B12 

A2 = (A22 - All) - A21 * B12 - (A33 - All) - (A44 - All) 

X + A43 * B34 

A3 = A(L1+1,L1) / B22 

GO TO 160 

C AD HOC SHIFT 

155 A1 = O.OEO 
A2 = 1.0E0 
A3 = 1.1605E0 
160 ITS = ITS + 1 
ITN = ITN - 1 
IF (.NOT. MATZ) LORI = LD 

C MAIN LOOP 

DO 260 K = L, NA 

NOTLAS = K .NE. NA .AND. ISH .EQ. 2 
K1 = K + 1 
K2 = K + 2 


EISP6932 

EISP693 3 

EISP6934 

EISP6935 

EISP6936 

EISP6937 

EISP6938 

EISP6939 

EISP6940 

EISP6941 

EISP6942 

EISP6943 

EISP6944 

EISP6945 

EISP6946 

EISP6947 

EISP6948 

EISP6949 

EISP6950 

EISP6951 

EISP6952 

EISP6953 

EISP6954 

EISP6955 

EISP6956 

EISP6957 

EISP6958 

EISP6959 

EISP6960 

EISP6961 

EISP6962 

EISP6963 

EISP6964 

EISP6965 

EISP6966 

EISP6967 

EISP6968 

EISP6969 

EISP697 0 

EISP6971 

EISP6972 

EISP6973 

EISP6974 

EISP6975 

EISP6976 

EISP6977 

EISP6978 

EISP6979 

EISP6980 

EISP6981 

EISP6982 

EISP6983 

EISP6984 

EISP6985 

EISP6986 

EISP6987 

EISP6988 

EISP6989 

EISP6990 

EISP6991 



-C 


170 


C 


180 

C 


— c 

190 


200 


210 

~C 


c 


KM1 = MAXO (K-l , L) 

LL = MINO (EN , Kl+ISH) 

IF (NOTLAS) GO TO 190 

ZERO A (K+l , K-l ) 

IF (K .EQ. L) GO TO 170 
A1 = A(K,KM1) 

A2 = A (K1 , KM1) 

S = ABS(Al) + ABS ( A2 ) 

IF (S .EQ. 0.0E0) GO TO 70 
U1 = A1 / S 

U2 = A2 / S 

R = SIGN (SQRT (U1*U1+U2*U2 ) ,U1) 
VI = -(Ul + R) / R 
V2 = -U2 / R 

U2 = V2 / VI 

DO 180 J = KM1 , ENORN 

T = A (K , J) + U2 * A (K1 , J) 

A (K, J) = A (K, J) + T * VI 
A (K1 , J) = A (K1 , J) + T * V2 
T = B (K, J) + U2 * B (K1 , J) 

B (K , J) = B (K, J) + T * VI 
B (K1 , J) = B (K1 , J) + T * V2 
CONTINUE 


IF (K .NE. L) A (K1 , KM1 ) = 0.0E0 
GO TO 240 

ZERO A(K+1,K-1) AND A(K+2,K-1) 

IF (K .EQ. L) GO TO 200 
A1 = A (K, KM1) 

A2 = A (K1 , KM1) 

A3 = A(K2,KM1) 

S = ABS(Al) + ABS (A2 ) + ABS (A3) 

IF (S .EQ. 0.0E0) GO TO 260 
Ul = A1 / S 
U2 = A2 / S 
U3 = A3 / S 

R = SIGN (SQRT (U1*U1+U2 *U2+U3 *U3 ) , Ul) 

VI = -(Ul + R) / R 

V2 = -U2 / R 

V3 = -U3 / R 

U2 = V2 / VI 

U3 = V3 / VI 


DO 210 J = KM1 , ENORN 

T = A (K , J) + U2 * A (K1 , J) + 
A (K , J) = A (K , J) + T * VI 
A (K1 , J) = A (K1 , J) + T * V2 

A (K2 , J) = A (K2 , J) + T * V3 

T = B (K, J) + U2 * B(K1,J) + 
B (K, J) = B (K, J) + T * VI 
B (K1 , J) = B (K1 , J) + T * V2 

B (K2 , J) = B (K2 , J) + T * V3 

CONTINUE 


U3 * A (K2 , J) 


U3 * B (K2 , J) 


IF (K .EQ. L) GO TO 220 
A (K1 , KM1) = 0.0E0 
A (K2 , KM1) = 0.0E0 

ZERO B (K+2 , K+l ) AND B(K+2,K) 

S = ABS (B (K2 , K2 ) ) + ABS (B (K2 , Kl) ) + ABS(B(K2,K)) 


EISP6992 

EISP6993 

EISP6994 

EISP6995 

EISP6996 

EISP6997 

EISP6998 

EISP6999 

EISP7000 

EISP7001 

EISP7002 

EISP7003 

EISP7004 

EISP7005 

EISP7006 

EISP7007 

EISP7008 

EISP7009 

EISP7010 

EISP7011 

EISP7012 

EISP7013 

EISP7014 

EISP7015 

EISP7016 

EISP7017 

EISP7018 

EISP7019 

EISP7020 

EISP7021 

EISP7022 

EISP7023 

EISP7024 

EISP7025 

EISP7026 

EISP7027 

EISP7028 

EISP7029 

EISP7030 

EISP703 1 

EISP7032 

EISP7033 

EISP7034 

EISP7035 

EISP7036 

EISP7037 

EISP7038 

EISP7039 

EISP7040 

EISP7041 

EISP7042 

EISP7043 

EISP7044 

EISP7045 

EISP7046 

EISP7047 

EISP7048 

EISP7049 

EISP7050 

EISP7051 


220 




IF (S .EQ. O.OEO) GO TO 240 

EISP7052 



Ul = B(K2,K2) / S 

EISP7053 



U2 = B(K2,K1) / S 

EISP7054 

— 


U3 = B(K2,K) / S 

EISP7055 



R = SIGN ( SQRT (U1*U1+U2*U2+U3*U3 ) , Ul) 

EISP7056 



VI = -(Ul + R) / R 

EISP7057 



V2 = -U2 / R 

EISP7058 



V3 = -U3 / R 

EISP7059 



U2 = V2 / VI 

EISP7060 



U3 = V3 / VI 

EISP7061 

-c 



EISP7062 



DO 230 I = LORI, LL 

EISP7063 



T = A ( I , K2 ) + U2 * A (I , Kl) + U3 * A(I,K) 

EISP7064 



A ( I , K2 ) = A ( I , K2 ) + T * VI 

EISP7065 



A (I , Kl) = A(I ,K1) + T * V2 

EISP7066 



A ( I , K) = A ( I , K) + T * V3 

EISP7067 



T = B (I , K2 ) + U2 * B (I , Kl) + U3 * B(I,K) 

EISP7068 



B ( I , K2 ) = B ( I , K2 ) + T * VI 

EISP7069 



B ( I , Kl) = B ( I , Kl ) + T * V2 

EISP7070 



B ( I , K) = B ( I , K) + T * V3 

EISP7071 

— 

230 

CONTINUE 

EISP7072 

c 



EISP7073 



B(K2,K) = O.OEO 

EISP7074 



B (K2 , Kl) = O.OEO 

EISP7075 



IF (.NOT. MATZ) GO TO 240 

EISP7076 

c 



EISP7077 



DO 235 I = 1, N 

EISP7078 

— 


T = Z ( I , K2 ) + U2 * Z (I ,K1) + U3 * Z(I,K) 

EISP7079 



Z (I ,K2) = Z(I,K2) + T * VI 

EISP7080 



Z (I , Kl) = Z ( I , Kl) + T * V2 

EISP7081 




Z (I , K) ■ Z ( I , K) + T * V3 

EISP7082 


235 

CONTINUE 

EISP7083 

c 

. 

ZERO B (K+l , K) 

EISP7084 


240 

S = ABS (B (Kl , Kl) ) + ABS (B (Kl , K) ) 

EISP7085 



IF (S .EQ. O.OEO) GO TO 260 

EISP7086 



Ul = B (Kl , Kl) / S 

EISP7087 



U2 = B(K1,K) / S 

EISP7088 

— 


R = SIGN (SQRT (U1*U1+U2*U2 ) ,U1) 

EISP7089 



VI = -(Ul + R) / R 

EISP7090 



V2 = -U2 / R 

EISP7091 




U2 = V2 / VI 

EISP7092 

c 



EISP7093 



DO 250 I = LORI, LL 

EISP7094 



T = A ( I , Kl) + U2 * A ( I , K) 

EISP7095 



A (I , Kl) = A(I ,K1) + T * VI 

EISP7096 



A (I , K) = A ( I , K) + T * V2 

EISP7097 



T = B(I ,K1) + U2 * B ( I , K) 

EISP7098 

— 


B (I , Kl ) = B (I , Kl) + T * VI 

EISP7099 



B (I , K) = B (I , K) + T * V2 

EISP7100 


250 

CONTINUE 

EISP7101 

n 





EISP7102 



B (Kl , K) = O.OEO 

EISP7103 



IF (.NOT. MATZ) GO TO 260 

EISP7104 

n 

'C. 



EISP7105 



DO 255 I = 1, N 

EISP7106 



T = Z (I , Kl) + U2 * Z ( I , K) 

EISP7107 



Z ( I , Kl) = Z ( I , Kl) + T * VI 

EISP7108 

— 


Z ( I , K) = Z ( I , K) + T * V2 

EISP7109 


255 

CONTINUE 

EISP7110 

r» 

Co 



EISP7111 



260 CONTINUE 


END QZ STEP 


C 

-C 

C 


GO TO 70 


SET ERROR — ALL EIGENVALUES HAVE NOT 
CONVERGED AFTER 30*N ITERATIONS 


1000 I ERR = EN 

C SAVE EPSB FOR USE BY QZVAL AND QZVEC 

1001 IF (N .GT. 1) B (N , 1) = EPSB 
RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 



ROUTINE NAME 
FROM EISPACK 


LATEST REVISION 


PURPOSE 


USAGE 

ARGUMENTS NM 


N 

A 


B 


PF262=QZVAL 


AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 


THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN QUASI -TRIANGULAR 
FORM AND THE OTHER IN UPPER TRIANGULAR FORM. 
IT REDUCES THE QUASI -TRIANGULAR MATRIX 
FURTHER, SO THAT ANY REMAINING 2 -BY- 2 BLOCKS 
CORRESPOND TO PAIRS OF COMPLEX EIGENVALUES, 
AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE 
GENERALIZED EIGENVALUES. IT IS USUALLY 
PRECEDED BY QZHES(PF260) AND QZIT(PF261) AND 
MAY BE FOLLOWED BY QZVEC (PF263 ) . 


CALL QZVAL (NM , N , A , B , ALFR, ALFI , BETA , MATZ , Z ) 

ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT. 

ON INPUT N IS THE ORDER OF THE MATRICES. 

ON INPUT A CONTAINS A REAL UPPER QUASI- 
TRIANGULAR MATRIX. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT A HAS BEEN REDUCED FURTHER TO A 
QUASI -TRIANGULAR MATRIX IN WHICH ALL NONZERO 
SUBDIAGONAL ELEMENTS CORRESPOND TO PAIRS OF 
COMPLEX EIGENVALUES. 

ON INPUT B CONTAINS A REAL UPPER TRIANGULAR 
MATRIX. 

MUST BE OF DIMENSION NM X N. 

IN ADDITION, LOCATION B(N,1) CONTAINS THE 
TOLERANCE QUANTITY (EPSB) COMPUTED AND SAVED 
IN QZIT (PF261) . 

ON OUTPUT B IS STILL IN UPPER TRIANGULAR 
FORM, ALTHOUGH ITS ELEMENTS HAVE BEEN ALTERED 
B (N, 1) IS UNALTERED. 


EISP7112 
EISP7113 
EISP7114 
EISP7115 
EISP7116 
EISP7117 
EISP7118 
EISP7119 
EISP7120 
EISP7121 
EISP7122 
QZVAL 2 
QZVAL 3 
QZVAL 4 
QZVAL 5 
QZVAL 6 
QZVAL 7 
QZVAL 8 
QZVAL 9 
QZVAL 10 
QZVAL 11 
QZVAL 12 
QZVAL 13 
QZVAL 14 
QZVAL 15 
QZVAL 16 
QZVAL 17 
QZVAL 18 
QZVAL 19 
QZVAL 20 
QZVAL 21 
QZVAL 22 
QZVAL 23 
QZVAL 24 
QZVAL 25 
QZVAL 26 
QZVAL 27 
QZVAL 28 
QZVAL 29 
QZVAL 30 
QZVAL 31 
QZVAL 32 
QZVAL 33 
QZVAL 34 
QZVAL 35 
QZVAL 36 
QZVAL 37 
QZVAL 38 
QZVAL 39 
QZVAL 40 
QZVAL 41 
QZVAL 42 
QZVAL 43 
QZVAL 44 
QZVAL 45 
QZVAL 46 
QZVAL 47 
QZVAL 48 
QZVAL 49 
QZVAL 50 



~c 

c 

c 

-c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 

_c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 

_c 

c 

c 

c 

~c 

c 

c 


c 

n 

n 

-c 

c 


c 

n 

-c 

c 


c 

r 


ALFR - ON OUTPUT ALFR CONTAINS THE REAL PART OF THE 
DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX 
THAT WOULD BE OBTAINED IF A WERE REDUCED 
COMPLETELY TO TRIANGULAR FORM BY UNITARY 
TRANSFORMATIONS. NON-ZERO VALUES OF ALFI 
OCCUR IN PAIRS, THE FIRST MEMBER POSITIVE AND 
THE SECOND NEGATIVE. 

MUST BE OF DIMENSION N. 

ALFI - ON OUTPUT ALFI CONTAINS THE IMAGINARY PART 


BETA 


MATZ 


REQUIRED ROUTINES 


QZVAL 51 
QZVAL 52 
QZVAL 53 
QZVAL 54 
QZVAL 55 
QZVAL 56 
QZVAL 57 
QZVAL 58 
QZVAL 59 
QZVAL 60 
QZVAL 61 
QZVAL 62 


OF THE DIAGONAL ELEMENTS OF OF THE TRIANGULAR QZVAL 63 


MATRIX THAT WOULD BE OBTAINED IF A WERE 
REDUCED COMPLETELY TO TRIANGULAR FORM BY 
UNITARY TRANSFORMATIONS. NON-ZERO VALUES 
OF ALFI OCCUR IN PAIRS, THE FIRST MEMBER 
POSITIVE AND THE SECOND NEGATIVE. 

MUST BE OF DIMENSION N. 


QZVAL 64 
QZVAL 65 
QZVAL 66 
QZVAL 67 
QZVAL 68 
QZVAL 69 
QZVAL 70 

ON OUTPUT BETA CONTAINS THE DIAGONAL ELEMENTS QZVAL 71 
OF THE CORRESPONDING B, NORMALIZED TO BE REAL QZVAL 72 
AND NON-NEGATIVE. THE GENERALIZED EIGENVALUESQZVAL 73 


ARE THEN THE RATIOS ( (ALFR+I*ALFI) /BETA) 
MUST BE OF DIMENSION N. 


ON INPUT MATZ SHOULD BE SET TO .TRUE. IF 
THE RIGHT HAND TRANSFORMATIONS ARE TO BE 
ACCUMULATED FOR LATER USE IN COMPUTING 
EIGENVECTORS, AND TO .FALSE. OTHERWISE. 


ON 

TO 


INPUT Z 
. TRUE . , 


CONTAINS, IF MATZ HAS BEEN SET 


QZVAL 74 
QZVAL 75 
QZVAL 76 
QZVAL 77 
QZVAL 78 
QZVAL 79 
QZVAL 80 
QZVAL 81 
QZVAL 82 
QZVAL 83 


IN THE REDUCTIONS BY QZHES(PF260) AND QZIT 
(PF261) IF PERFORMED, OR ELSE THE IDENTITY 
MATRIX. IF MATZ HAS BEEN SET TO .FALSE., Z 
IS NOT REFERENCED. 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT Z CONTAINS THE PRODUCT OF THE 
RIGHT HAND TRANSFORMATIONS (FOR ALL THREE 
STEPS) IF MATZ HAS BEEN SET TO .TRUE. 

NONE 


REMARKS 1. 


EXAMPLE : 

PROGRAM TQZVAL (OUTPUT , TAPE6=OUTPUT) 

DIMENSION A(5,5) ,B(5,5) , ALFR (5) , ALFI (5) , BETA (5) ,Z(5,5) 
LOGICAL MATZ 


N = 5 
NM = 5 
MATZ = 
EPS1 = 


. TRUE . 
0.0E0 


THE TRANSFORMATION MATRIX PRODUCED QZVAL 84 

QZVAL 85 
QZVAL 86 
QZVAL 87 
QZVAL 88 
QZVAL 89 
QZVAL 90 
QZVAL 91 
QZVAL 92 
QZVAL 93 
QZVAL 94 
QZVAL 95 
QZVAL 96 
QZVAL 97 
QZVAL 98 
QZVAL 99 
QZVAL100 
QZVAL101 
QZVAL102 
QZVAL103 
QZVAL104 
QZVAL105 
QZVAL106 
QZVAL107 
QZVAL108 
QZVAL109 
QZVAL110 


THIS SUBROUTINE IS THE THIRD STEP OF THE 
ALGORITHM FOR SOLVING GENERALIZED MATRIX 
EIGENVALUE PROBLEMS, SIAM J. NUMER. ANAL. 
241-256(1973) BY MOLER AND STEWART. 


QZ 


10 , 



QZVAL111 


c 

DATA A /10 

• f 2 • f 3 • , 2*1 . , 2 . , 12 . f 1 • f 2 • f 1 • f 3 • f 1 • / 1 1 • / 

QZVAL112 

c 

* 1. , 

-1. ,1. ,2. , 

1. ,9. ,3*1. ,-l. ,1. ,15. / 

QZVAL113 

-c 



QZVAL114 

c 

DATA B / 12 

. ,1. ,-l. ,2 

. ,2*1. ,14. ,1. , — 1 . , 1 . ,— 1. ,1* , 

QZVAL115 

c 

* 16. 

,-l. ,1. ,2. 

,-l. ,-l. ,12. ,-l. ,3*1. ,-l. ,11. / 

QZVAL116 

c 




QZVAL117 

c 

CALL QZHES (NM, N, A, B 

,MATZ, Z) 

QZVAL118 

c 

CALL QZIT(NM,N, A, B, 

EPS1 , MATZ , Z , IERR) 

QZVAL119 

c 

CALL QZVAL(NM,N, A, B 

, ALFR , ALFI , BETA , MATZ , Z ) 

QZVAL120 

-C 

WRITE (6, 99) I ERR 


QZVAL121 

c 

WRITE ( 6 , 100) ALFR, ALFI, BETA, ( (Z(I,J) ,1=1,5) ,J=1,5) 

QZVAL122 

C99 

FORMAT ( 1H1 

,8H I ERR = 

,14) 

QZVAL123 

_C100 

FORMAT ( 1H 

,8H ALFR = 

/1H , 5 (G8 . 2 , 2X) / 

QZVAL124 

C 

* 8H 

ALFI = /1H 

,5(G8.2,2X) / 

QZVAL125 

C 

* 8H 

BETA = /1H 

, 5 (G8 . 2 , 2X) / 

QZVAL126 

c 

* 5H 

Z = / 5 ( 1H 

, 5 (G8 . 2 , 2X) / ) ) 

QZVAL127 

c 

STOP 



QZVAL128 

c 

END 



QZVAL129 

c 




QZVAL130 

-c 

OUTPUT : 



QZVAL131 

c 




QZVAL132 

c 

I ERR = 

0 


QZVAL133 

_c 

ALFR = 



QZVAL134 

c 

15. 

7.2 

16. 10. 8.6 

QZVAL135 

c 

ALFI = 



QZVAL136 

c 

0. 

0. 

• 

o 

• 

o 

o 

QZVAL137 

~C 

BETA = 



QZVAL138 

c 

9.9 

17. 

14. 11. 13. 

QZVAL139 

c 

Z = 



QZVAL140 

-C 

.24 

-. 54E-01 

.21 -.27 -.91 

QZVAL141 

c 

-.54 

.25 

.65 -.46 .13 

QZVAL142 

c 

.49 

.56 

.49 .45 . 75E-01 

QZVAL143 

c 

-.60 

.48 

-.29 .44 -.38 

QZVAL144 

c 

-.25 

-.63 

.45 .57 -. 94E-01 

QZVAL145 


C 

C- 


SUBROUTINE Q Z VAL ( NM , N , A , B , ALFR , ALFI , BETA , MATZ , Z ) 


QZVAL146 

QZVAL147 

EISP7123 


implicit real*8 (a-h,o-z) 

_ INTEGER I, J,N, EN,NA,NM,NN, ISW 

REAL* 8 A (NM, N) , B (NM, N) , ALFR(N) , ALFI (N) , BETA(N) ,Z(NM,N) 

REAL* 8 C, D, E , R, S , T, AN, A1 , A2 , BN, CQ, CZ , DI , DR, El , TI , TR, U1 , 

X U2 , VI , V2 , All , All , A12 , A2I , A21 , A22 , Bll , B12 , B22 , SQI , SQR, 

X SSI , SSR, SZI , SZR, A11I , A11R, A12I , A12R, A22I , A22R, EPSB 

LOGICAL MATZ 
EPSB = B (N, 1) 

- ISW = 1 

C FIND EIGENVALUES OF QUASI -TRIANGULAR MATRICES. 

C FOR EN=N STEP -1 UNTIL 1 DO — 

_ DO 510 NN = 1, N 

EN = N + 1 - NN 
NA = EN - 1 

IF (ISW .EQ. 2) GO TO 505 
IF (EN .EQ. 1) GO TO 410 
IF (A (EN, NA) .NE. 0.0E0) GO TO 420 
C 1-BY-l BLOCK, ONE REAL ROOT 

- 410 ALFR(EN) = A(EN, EN) 

IF (B (EN, EN) .LT. 0.0E0) ALFR(EN) = -ALFR(EN) 

BETA ( EN ) = ABS (B (EN , EN) ) 


EISP7124 

EISP7125 

EISP7126 

EISP7127 

EISP7128 

EISP7129 

EISP7130 

EISP7131 

EISP7132 

EISP7133 

EISP7134 

EISP7135 

EISP7136 

EISP7137 

EISP7138 

EISP7139 

EISP7140 

EISP7141 

EISP7142 

EISP7143 

EISP7144 



o o 


c 


c 


_c 


ALFI(EN) = O.OEO 
GO TO 510 

2-BY-2 BLOCK 

420 IF (ABS (B (NA,NA) ) .LE. EPSB) GO TO 455 
IF (ABS (B (EN, EN) ) .GT. EPSB) GO TO 430 
A1 = A ( EN , EN ) 

A2 = A(EN,NA) 

BN = O.OEO 
GO TO 435 

430 AN = ABS (A (NA, NA) ) + ABS (A (NA, EN) ) + ABS (A (EN, NA) ) 

X + ABS (A(EN, EN) ) 

BN = ABS (B (NA, NA) ) + ABS (B (NA, EN) ) + ABS (B (EN, EN) ) 

All = A (NA, NA) / AN 

A12 = A(NA, EN) / AN 

A21 = A(EN,NA) / AN 

A22 = A ( EN , EN ) / AN 

Bll = B (NA, NA) / BN 

B12 = B (NA, EN) / BN 

B22 = B (EN, EN) / BN 

E = All / Bll 

El = A22 / B22 

S = A21 / (Bll * B22 ) 

T = (A22 - E * B22) / B22 

IF (ABS (E) .LE. ABS (El)) GO TO 431 

E = El 

T = (All - E * Bll) / Bll 

431 C = 0.5E0 * (T - S * B12) 

D=C*C+S* (A12 - E * B12 ) 

IF (D .LT. O.OEO) GO TO 480 

TWO REAL ROOTS. 

ZERO BOTH A(EN,NA) AND B (EN,NA) 

E = E + (C + SIGN ( SQRT (D) , C) ) 


All = All 

- E 

* 

Bll 

A12 = A12 

- E 

* 

B12 

A22 = A22 

- E 

* 

B22 


IF (ABS (All) + ABS (A12 ) .LT. 

X ABS (A21) + ABS (A22 ) ) GO TO 432 

A1 = A12 
A2 = All 
GO TO 435 
432 A1 = A22 
A2 = A21 

CHOOSE AND APPLY REAL Z . . 

435 S = ABS(Al) + ABS(A2) 

U1 = A1 / S 
U2 = A2 / S 

R = SIGN (SQRT (U1*U1+U2*U2 ) ,U1) 

VI = -(Ul + R) / R 
V2 = -U2 / R 
U2 = V2 / VI 

DO 440 1=1, EN 

T = A ( I , EN) + U2 * A ( I , NA) 

A (I , EN) = A ( I , EN ) + T * VI 

A ( I , NA) = A (I , NA) + T * V2 

T = B(I,EN) + U2 * B ( I , NA) 

B ( I , EN ) = B ( I , EN) + T * VI 

B ( I , NA) = B ( I , NA) + T * V2 

440 CONTINUE 
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EISP7201 

EISP7202 

EISP7203 

EISP7204 


C 



— 


IF (.NOT. MATZ) GO TO 450 

EISP7205 

c 



EISP7206 



DO 445 I = 1, N 

EISP7207 

— 


T = Z ( I , EN ) + U2 * Z (I , NA) 

EISP7208 



Z ( I , EN ) = Z ( I , EN ) + T * VI 

EISP7209 



Z ( I , NA) = Z ( I , NA) + T * V2 

EISP7210 


445 

CONTINUE 

EISP7211 

~c 



EISP72 12 


450 

IF (BN .EQ. 0.0E0) GO TO 475 

EISP7213 



IF (AN .LT. ABS(E) * BN) GO TO 455 

EISP7214 

— 


A1 = B (NA, NA) 

EISP7215 



A2 = B (EN, NA) 

EISP7216 



GO TO 460 

EISP7217 


455 

A1 = A (NA, NA) 

EISP7218 



A2 = A ( EN , NA ) 

EISP7219 

c 

# f 

CHOOSE AND APPLY REAL Q 

EISP7220 


460 

S = ABS(Al) + ABS ( A2 ) 

EISP7221 



IF (S .EQ. 0.0E0) GO TO 475 

EISP7222 



U1 = A1 / S 

EISP7223 



U2 = A2 / S 

EISP7224 

— 


R = SIGN ( SQRT (U1*U1+U2 *U2 ) ,U1) 

EISP7225 



VI = -(Ul + R) / R 

EISP7226 



V2 = -U2 / R 

EISP7227 



U2 = V2 / VI 

EISP7228 

c 



EISP7229 



DO 470 J = NA, N 

EISP7230 



T = A (NA, J) + U2 * A ( EN , J ) 

EISP7231 



A (NA, J) = A (NA, J) + T * VI 

EISP7232 



A ( EN , J ) = A ( EN , J ) + T * V2 

EISP723 3 



T = B (NA, J) + U2 * B(EN, J) 

EISP7234 

— 


B (NA, J) = B (NA, J) + T * VI 

EISP7235 



B ( EN , J ) = B ( EN , J ) + T * V2 

EISP7236 


470 

CONTINUE 

EISP7237 

_c 



EISP7238 


475 

A(EN,NA) = 0.0E0 

EISP7239 



B (EN, NA) = 0.0E0 

EISP7240 



ALFR(NA) = A (NA, NA) 

EISP7241 

— 


ALFR(EN) = A ( EN , EN ) 

EISP7242 



IF (B (NA, NA) .LT. 0.0E0) ALFR(NA) = -ALFR(NA) 

EISP7243 



IF ( B ( EN , EN ) .LT. 0.0E0) ALFR(EN) = -ALFR(EN) 

EISP7244 




BETA ( NA ) = ABS (B (NA, NA) ) 

EISP7245 



BETA(EN) = ABS (B (EN, EN) ) 

EISP7246 



ALFI (EN) = 0.0E0 

EISP7247 



ALFI (NA) = 0.0E0 

EISP7248 



GO TO 505 

EISP7249 

c 

. , 

TWO COMPLEX ROOTS 

EISP7250 


480 

E = E + C 

EISP7251 

— 


El = SQRT(-D) 

EISP7252 



A11R = All - E * Bll 

EISP7253 



A11I = El * Bll 

EISP7254 



A12R = A12 - E * B12 

EISP7255 



A12I = El * B12 

EISP7256 



A22R = A22 - E * B22 

EISP7257 



A22I = El * B22 

EISP7258 



IF (ABS(AllR) + ABS(AllI) + ABS(A12R) + ABS(A12I) .LT. 

EISP7259 


X 

ABS (A21) + ABS (A22R) + ABS(A22I)) GO TO 482 

EISP7260 



A1 = A12R 

EISP7261 

— 


All = A12I 

EISP7262 



A2 = -A11R 

EISP7263 



A2I = -A11I 

EISP7264 



482 


GO TO 485 


A1 = A22R 
All = A22I 
A2 = -A2 1 
A2I = 0.0E0 

C CHOOSE COMPLEX Z 

485 CZ = SQRT (A1*A1+A1I*A1I) 

IF (CZ .EQ. 0.0E0) GO TO 487 
SZR = (A1 * A2 + All * A2I ) / CZ 

SZI = (A1 * A2I - All * A2) / CZ 

R = SQRT (CZ*CZ+SZR*SZR+SZI*SZI) 
CZ = CZ / R 
SZR = SZR / R 

SZI = SZI / R 


GO TO 490 

487 SZR = l.OEO 
SZI = O.OEO 

490 IF (AN .LT. (ABS(E) + El) * BN) GO TO 492 
A1 = CZ * Bll + SZR * B12 
All = SZI * B12 
A2 = SZR * B22 
A2I = SZI * B22 
GO TO 495 


__ 492 A1 = CZ * All + SZR * A12 

All = SZI * A12 
A2 = CZ * A21 + SZR * A22 
A2I = SZI * A22 

~C CHOOSE COMPLEX Q 

495 CQ = SQRT (A1*A1+A1I*A1I ) 

IF (CQ .EQ. O.OEO) GO TO 497 
- SQR = (A1 * A2 + All * A2I) / CQ 

SQI = (A1 * A2I - All * A2) / CQ 
R = SQRT (CQ*CQ+SQR*SQR+SQI*SQI) 


CQ = CQ / R 
SQR = SQR / R 
SQI = SQI / R 
GO TO 500 

- 497 SQR = l.OEO 

SQI = O.OEO 

C COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT 

_C IF TRANSFORMATIONS WERE APPLIED 

500 SSR = SQR * SZR + SQI * SZI 
SSI = SQR * SZI - SQI * SZR 
1 = 1 

TR = CQ * CZ * All + CQ * SZR * A12 + SQR * CZ * A21 
X + SSR * A22 

TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22 

- DR = CQ * CZ * Bll + CQ * SZR * B12 + SSR * B22 
DI = CQ * SZI * B12 + SSI * B22 

GO TO 503 
_ 502 1=2 

TR = SSR * All - SQR * CZ * A12 - CQ * SZR * A21 
X + CQ * CZ * A22 

TI = -SSI * All - SQI * CZ * A12 + CQ * SZI * A21 
DR = SSR * Bll - SQR * CZ * B12 + CQ * CZ * B22 
DI = -SSI * Bll - SQI * CZ * B12 
503 T = TI * DR - TR * DI 
J = NA 

IF (T .LT. O.OEO) J = EN 
R = SQRT ( DR*DR+DI *DI ) 
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EISP7322 
EISP7323 
EISP7324 



BETA ( J ) = BN 
ALFR(J) = AN 
ALFI(J) = AN 
IF (I .EQ. 1) 


* R 

* (TR * 

* T / R 
GO TO 502 


DR + TI * DI) / R 


505 

510 


ISW = 
CONTINUE 
B (N f 1) = 


3 - ISW 


EPSB 


C** 

C 

C 

_C 

C 

C 

_C 

C 

C 

C 

-C 

C 

C 

c 
c 
c 
~c 
c 
c 
-c 
c 
c 
_c 
c 
c 
c 
— c 
c 
c 
_c 
c 
c 
c 
c 
c 
c 
-c 
c 
c 

c 

c 

_c 

c 

c 

c 

-c 

c 

c 


RETURN 

THIS PROGRAM VALID 
END 

ROUTINE NAME 
FROM EISPACK 


ON FTN4 AND FTN5 ** 
- PF2 63=QZVEC 


LATEST REVISION 


PURPOSE 


USAGE 

ARGUMENTS NM 


N 


B 


ALFR 


AUGUST 1,1984 
COMPUTER SCIENCES CORP. 


HAMPTON, VA. 


- THIS SUBROUTINE ACCEPTS A PAIR OF REAL 
MATRICES, ONE OF THEM IN QUASI -TRIANGULAR 
FORM (IN WHICH EACH 2-BY-2 BLOCK CORRESPONDS 
TO A PAIR OF COMPLEX EIGENVALUES) AND THE 
OTHER IN UPPER TRIANGULAR FORM. IT COMPUTES 
THE EIGENVECTORS OF THE TRIANGULAR PROBLEM 
AND TRANSFORMS THE RESULTS BACK TO THE 
ORIGINAL COORDINATE SYSTEM. IT IS USUALLY 
PRECEDED BY QZHES (PF260) , QZIT(PF261), AND 
QZVAL ( PF2 62) . 

- CALL QZVEC (NM,N, A, B, ALFR, ALFI , BETA, Z) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT . 

- ON INPUT N IS THE ORDER OF THE MATRICES. 

- ON INPUT A CONTAINS A REAL UPPER QUASI- 
TRI ANGULAR MATRIX. 

MUST BE OF DIMENSION NM X N . 


EISP7325 
EISP7326 
EISP7327 
EISP7328 
EISP7329 
EISP7330 
EISP7331 
EISP7332 
EISP7333 
EISP7334 
EISP7335 
QZVEC 2 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 21 
QZVEC 22 
QZVEC 23 
QZVEC 24 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 30 
QZVEC 31 
QZVEC 32 
QZVEC 
QZVEC 
QZVEC 
QZVEC 


3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 


25 

26 

27 

28 
29 


ON OUTPUT A IS UNALTERED. ITS SUBDIAGONAL 
ELEMENTS PROVIDE INFORMATION ABOUT THE STORAGEQZVEC 

OF THE COMPLEX EIGENVECTORS. QZVEC 

QZVEC 

ON INPUT B CONTAINS A REAL UPPER TRIANGULAR QZVEC 


33 

34 

35 

36 

37 

38 

39 

40 


MATRIX. IN ADDITION, LOCATION B(N,1) CONTAINSQZVEC 41 

w m i. t ^ fT X TT? O ^ O 


nnxx\j-/\t * ... 

THE TOLERANCE QUANTITY (EPSB) COMPUTED AND 
SAVED IN QZIT (PF261) . 

MUST BE OF DIMENSION NM X N. 

ON OUTPUT B HAS BEEN DESTROYED. 

ON INPUT ALFR IS A VECTOR SUCH THAT THE 
RATIOS ( (ALFR+I*ALFI) /BETA) ARE THE 
GENERALIZED EIGENVALUES. THEY ARE USUALLY 


QZVEC 42 
QZVEC 43 
QZVEC 44 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 
QZVEC 


45 

46 

47 

48 

49 

50 


non 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

"C 

c 

c 

-C 

c 

c 

c 

'c 

c 

c 

-c 

c 

c 

_c 
c 
c 
c 
c 
c 
c 
. c 
c 
c 
_c 
c 
c 
c 
~c 
c 
c 
~-c 
c 
c 
c 


— c 
c 
c 
_c 
c 
c 
_c 
c 
c 
c 
-c 
c 
c 


ALFI 


BETA 


OBTAINED FROM QZVAL(PF262) . 
MUST BE OF DIMENSION N. 


QZVEC 51 
QZVEC 52 
QZVEC 53 

- ON INPUT ALFI IS A VECTOR SUCH THAT THE RATIOSQZVEC 54 

( (ALFR+I*ALFI) /BETA) ARE THE GENERALIZED QZVEC 55 

EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVEC 56 

QZVAL(PF262) . §zVEC58 

MUST BE OF DIMENSION N. QZVEC 59 

- ON INPUT BETA IS A VECTOR SUCH THAT THE RATIOSQZVEC 60 

( (ALFR+I*ALFI) /BETA) ARE THE GENERALIZED QZVEC 61 

EIGENVALUES. THEY ARE USUALLY OBTAINED FROM QZVEC 62 

QZVAL(PF262) . OZVEC64 

MUST BE OF DIMENSION N. Q2VEC 64 

QZVEC 65 

- ON INPUT Z CONTAINS THE TRANSFORMATION MATRIX QZVEC 66 
PRODUCED IN THE REDUCTIONS BY QZHES (PF260) , QZVEC 67 
QZIT (PF261) , AND QZVAL(PF262) , IF PERFORMED. QZVEC 68 
IF THE EIGENVECTORS OF THE TRIANGULAR PROBLEM QZVEC 69 


Z MUST CONTAIN THE IDENTITY 


ARE DESIRED, 

MATRIX. 

MUST BE OF DIMENSION NM X N. 


QZVEC 70 
QZVEC 71 
QZVEC 72 
QZVEC 73 
QZVEC 74 
QZVEC 75 


ON OUTPUT Z CONTAINS THE REAL AND IMAGINARY 
PARTS OF THE EIGENVECTORS. IF ALFI (I) .EQ. 

0.0, THE I-TH EIGENVALUE IS REAL AND THE I-TH QZVEC 76 
COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF QZVEC 77 

ALFI (I) .NE. 0.0, THE I-TH EIGENVALUE IS QZVEC 78 

COMPLEX. IF ALFI (I) -GT. 0.0, THE EIGENVALUE QZVEC 79 
IS THE FIRST OF A COMPLEX PAIR AND THE I-TH QZVEC 80 
AND (I+1)-TH COLUMNS OF Z CONTAIN ITS EIGEN- QZVEC 81 
VECTOR. IF ALFI ( I ) .LT. 0.0, THE EIGEN- QZVEC ?? 

VALUE IS THE SECOND OF A COMPLEX PAIR AND THEQZVEC ?? 


(I-l)-TH AND I-TH COLUMNS OF Z CONTAIN THE 
CONJUGATE OF ITS EIGENVECTOR. EACH EIGEN- 
VECTOR IS NORMALIZED SO THAT THE MODULUS 
OF ITS LARGEST COMPONENT IS 1.0 . 


REQUIRED ROUTINES - NONE 


REMARKS 1. 


THIS SUBROUTINE IS THE OPTIONAL FOURTH STEP 
OF THE QZ ALGORITHM FOR SOLVING GENERALIZED 
MATRIX EIGENVALUE PROBLEMS, SIAM J. NUMER. 
ANAL. 10, 241-256(1973) BY MOLER AND STEWART. 


EXAMPLE : 

PROGRAM TQZVEC (OUTPUT, TAPE6=OUTPUT) 

DIMENSION A ( 5 , 5) ,B(5,5) , ALFR ( 5 ) , ALFI (5) , BETA (5) ,Z(5,5) 
LOGICAL MATZ 


N = 5 
NM = 5 
MATZ = 
EPS1 = 


. TRUE . 
0.0E0 


DATA A / 10 . , 2 . , 3 . , 2*1 . , 2 . , 12 . , 1 . , 2 . , 1 . , 3 . , 1 
fc 1.,— l.,l.,2.,l.,9.,3*l.,“ l.,l.,15. 


, 11 ., 

/ 


DATA B /12 . , 1 . ,-l. , 2 . , 2*1. , 14 . , 1. , -1. , 1. ,-l- , I* , 


QZVEC 84 
QZVEC 85 
QZVEC 86 
QZVEC 87 
QZVEC 88 
QZVEC 89 
QZVEC 90 
QZVEC 91 
QZVEC 92 
QZVEC 93 
QZVEC 94 
QZVEC 95 
QZVEC 96 
QZVEC 97 
QZVEC 98 
QZVEC 99 
QZVEC100 
QZVEC101 
QZVEC102 
QZVEC103 
QZVEC104 
QZVEC105 
QZVEC106 
QZVEC107 
QZVEC108 
QZVEC109 
QZVEC110 


c 


* 16. ,-1. ,1. ,2. ,-l. ,-l. ,12. ,-l. ,3*1. ,-l. ,11. / 

QZVEC111 

c 



QZVEC112 

c 


CALL QZHES (NM, N, A, B,MATZ , Z) 

QZVEC113 

-c 


CALL QZIT (NM, N, A, B, EPS1 ,MATZ , Z , I ERR) 

QZVEC114 

c 


CALL QZVAL ( NM , N , A , B , ALFR , ALFI , BETA , MATZ , Z ) 

QZVEC115 

c 


CALL QZVEC (NM,N, A, B, ALFR, ALFI, BETA, Z) 

QZVEC116 

c 


WRITE (6, 99) I ERR 

QZVEC117 

c 


WRITE (6, 100) ( (Z(I,J) ,1=1,5) ,J=1,5) 

QZVEC118 

C99 

FORMAT (1H1,7HIERR = ,14) 

QZVEC119 

C100 

FORMAT ( 5H Z = /5(1H , 5 (G8 . 2 , 2X) / ) ) 

QZVEC120 

~c 


STOP 

QZVEC121 

c 


END 

QZVEC122 

c 



QZVEC123 

_c 


OUTPUT : 

QZVEC124 

c 



QZVEC125 

c 


I ERR = 0 

QZVEC126 

c 


Z = 

QZVEC127 

"c 


.26 - . 59E-01 .23 -.30 -1.0 

QZVEC128 

c 


-.85 .39 1.0 -.69 .26 

QZVEC129 

c 


1.0 1.0 .85 .88 . 54E-01 

QZVEC130 

-c 


-1.0 .83 -.39 .72 -.46 

QZVEC131 

c 


-.45 -.84 .65 1.0 - . 19E-01 

QZVEC132 

c 



QZVEC133 





c- 






SUBROUTINE QZVEC (NM , N , A , B , ALFR , ALFI , BETA , Z ) 

EISP7336 

U 


implicit real*8 (a-h,o-z) 

EISP7337 



INTEGER I , J , K , M , N , EN , I I , J J , NA , NM , NN , I SW , ENM2 

EISP7338 



REAL* 8 A (NM, N) , B ( NM , N ) , ALFR(N) , ALFI (N) , BETA(N) ,Z(NM,N) 

EISP7339 



REAL* 8 D,Q,R,S,T,W,X,Y, DI , DR, RA, RR, SA, TI , TR, T1 , T2 , W1 , XI , 

EISP7340 



X Z Z , Z 1 , ALFM , ALMI , ALMR , BETM , EPSB 

EISP7341 



EPSB = B (N, 1) 

EISP7342 



ISW = 1 

EISP7343 

c 


FOR EN-N STEP -1 UNTIL 1 DO — 

EISP7344 



DO 800 NN = 1, N 

EISP7345 



EN = N + 1 - NN 

EISP7346 



NA = EN - 1 

EISP7347 

— 


IF (ISW .EQ. 2) GO TO 795 

EISP7348 



IF (ALFI (EN) .NE. 0.0E0) GO TO 710 

EISP7349 

c 


REAL VECTOR 

EISP7350 



M = EN 

EISP7351 



B ( EN , EN ) = 1.0E0 

EISP7352 



IF (NA .EQ. 0) GO TO 800 

EISP7353 



ALFM = ALFR(M) 

EISP7354 



BETM = BETA (M) 

EISP7355 

c 


FOR I-EN-1 STEP -1 UNTIL 1 DO — 

EISP7356 



DO 700 II = 1, NA 

EISP7357 

— 


I = EN - II 

EISP7358 



W = BETM * A(I, I) - ALFM * B(I,I) 

EISP7359 



R = 0.0E0 

EISP7360 

_c 



EISP7361 



DO 610 J = M, EN 

EISP7362 


610 

R = R + (BETM * A ( I , J) - ALFM * B(I,J)) * B(J,EN) 

EISP7363 

c 



EISP7364 



IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 630 

EISP7365 



IF (BETM * A ( I , 1-1) .EQ. 0.0E0) GO TO 630 

EISP7366 



ZZ = W 

EISP7367 

-- 


S = R 

EISP7368 



GO TO 690 

EISP7369 


630 

M = I 

EISP7370 



IF (ISW .EQ. 2) GO TO 640 


C REAL 1-BY-l BLOCK 

T = W 

IF (W .EQ. 0.0E0) T = EPSB 
B ( I , EN ) = -R / T 
GO TO 700 

C REAL 2-BY-2 BLOCK 

640 X = BETM * A(I,I+1) - ALFM * B (1,1+1) 


Y = BETM * A (1+1,1) 

Q = W * ZZ — X * Y 
T = (X * S - ZZ * R) / Q 
B ( I , EN ) = T 

IF (ABS(X) .LE. ABS(ZZ)) GO TO 650 
B (1+1 , EN) = (-R - W * T) /X 
GO TO 690 


650 B ( 1+1 , EN) = (-S - Y * T) / ZZ 

690 ISW = 3 - ISW 

700 CONTINUE 

C END REAL VECTOR 

GO TO 800 

C COMPLEX VECTOR 

710 M = NA 


ALMR = ALFR(M) 
ALMI = ALFI(M) 
BETM = BETA(M) 


c LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT 

C EIGENVECTOR MATRIX IS TRIANGULAR 

Y = BETM * A (EN, NA) 


B (NA, NA) = -ALMI * B ( EN , EN ) / Y 

B (NA, EN) = (ALMR * B ( EN , EN ) - BETM * A(EN,EN) ) / Y 

B(EN,NA) = 0.0E0 

B (EN, EN) = 1.0E0 

ENM2 = NA - 1 

IF (ENM2 .EQ. 0) GO TO 795 

C FOR I=EN-2 STEP -1 UNTIL 1 DO — 

DO 790 II = 1, ENM2 
I = NA - II 

W = BETM * A (I, I) - ALMR * B(I,I) 

W1 = -ALMI * B (I , I) 

RA = 0.0E0 
SA = 0.0E0 
C 

DO 760 J = M, EN 

X = BETM * A(I,J) - ALMR * B(I,J) 

XI = -ALMI * B ( I , J) 

RA = RA + X * B ( J , NA) - XI * B(J,EN) 

SA = SA + X * B ( J, EN) + XI * B(J,NA) 

- 760 CONTINUE 

C 

IF (I .EQ. 1 .OR. ISW .EQ. 2) GO TO 770 
IF (BETM * A(I , 1-1) .EQ. 0.0E0) GO TO 770 
ZZ = W 
Z1 = W1 
R = RA 
S = SA 
ISW = 2 
GO TO 790 

- 770 M = I 

IF (ISW .EQ. 2) GO TO 780 

COMPLEX 1-BY-l BLOCK 


EISP7371 

EISP7372 

EISP7373 

EISP7374 

EISP7375 

EISP7376 

EISP7377 

EISP7378 

EISP7379 

EISP7380 

EISP7381 

EISP7382 

EISP7383 

EISP7384 

EISP7385 

EISP7386 

EISP7387 

EISP7388 

EISP7389 

EISP7390 

EISP7391 

EISP7392 

EISP7393 

EISP7394 

EISP7395 

EISP7396 

EISP7397 

EISP7398 

EISP7399 

EISP7400 

EISP7401 

EISP7402 

EISP7403 

EISP7404 

EISP7405 

EISP7406 

EISP7407 

EISP7408 

EISP7409 

EISP7410 

EISP7411 

EISP7412 

EISP7413 

EISP7414 

EISP7415 

EISP7416 

EISP74 17 

EISP7418 

EISP7419 

EISP7420 

EISP7421 

EISP7422 

EISP7423 

EISP7424 

EISP7425 

EISP7426 

EISP7427 

EISP7428 

EISP7429 

EISP7430 


C 


o on o o o ooo 


773 


TR = -RA 
TI = -SA 
DR = W 
DI = W1 


C COMPLEX DIVIDE (T1,T2) = (TR,TI) / (DR, DI) 

775 IF (ABS(DI) .GT. ABS(DR)) GO TO 777 


RR = DI / DR 
D = DR + DI * RR 
Tl = (TR + TI * RR) / D 

T2 = (TI - TR * RR) / D 

GO TO (787,782) , ISW 
CALL GOTOER 
777 RR = DR / DI 

D = DR * RR + DI 
Tl = (TR * RR + TI) / D 

T2 = (TI * RR - TR) / D 

GO TO (787,782) , ISW 
CALL GOTOER 


C COMPLEX 2 -BY- 2 BLOCK 

780 X = BETM * A(I,I+1) - ALMR * B(I,I+1) 

XI = -ALMI * B( 1,1+1) 

Y = BETM * A (1+1,1) 


TR = Y * RA - W * R + W1 * S 

TI = Y * SA - W * S - W1 * R 

DR = W * ZZ - W1 * Z1 - X * Y 
DI = W * Z1 + W1 * ZZ - XI * Y 

IF (DR .EQ. 0.0E0 .AND. DI .EQ. 0.0E0) DR = EPSB 
GO TO 775 

782 B ( 1+1 , NA) = Tl 

B ( 1+1 , EN) = T2 
ISW = 1 

IF (ABS(Y) .GT. ABS(W) + ABS(Wl)) GO TO 785 
TR = -RA - X * B (1+1 , NA) + XI * B(I+1,EN) 

TI = -SA - X * B (1+1 , EN) - XI * B(I+1,NA) 

GO TO 773 


785 Tl = (-R - ZZ * B(I+1,NA) + Z1 * B(I+1,EN)) / Y 

T2 = (-S - ZZ * B (1+1 , EN) - Z1 * B(I+1,NA)) / Y 
- 787 B ( I , NA) = Tl 

B ( I , EN ) = T2 
790 CONTINUE 

_C END COMPLEX VECTOR 

795 ISW = 3 - ISW 


800 CONTINUE 

END BACK SUBSTITUTION. 

TRANSFORM TO ORIGINAL COORDINATE SYSTEM. 

FOR J=N STEP -1 UNTIL 1 DO — 

DO 880 JJ = 1, N 
J = N + 1 - JJ 

DO 880 I = 1, N 
ZZ = 0.0E0 

DO 860 K = 1, J 

860 ZZ = ZZ + Z ( I , K) * B (K, J) 

Z(I,J) = ZZ 
880 CONTINUE 

NORMALIZE SO THAT MODULUS OF LARGEST 

COMPONENT OF EACH VECTOR IS 1. 

(ISW IS 1 INITIALLY FROM BEFORE) 


EISP7431 

EISP7432 

EISP7433 

EISP7434 

EISP7435 

EISP7436 

EISP7437 

EISP7438 

EISP7439 

EISP7440 

EISP7441 

EISP7442 

EISP7443 

EISP7444 

EISP7445 

EISP7446 

EISP7447 

EISP7448 

EISP7449 

EISP7450 

EISP7451 

EISP7452 

EISP7453 

EISP7454 

EISP7455 

EISP7456 

EISP7457 

EISP7458 

EISP7459 

EISP7460 

EISP7461 

EISP7462 

EISP7463 

EISP7464 

EISP7465 

EISP7466 

EISP7467 

EISP7468 

EISP7469 

EISP7470 

EISP7471 

EISP7472 

EISP7473 

EISP7474 

EISP7475 

EISP7476 

EISP7477 

EISP7478 

EISP7479 

EISP7480 

EISP7481 

EISP7482 

EISP7483 

EISP7484 

EISP7485 

EISP7486 

EISP7487 

EISP7488 

EISP7489 

EISP7490 



o o 


890 


C 

_C 


DO 950 J = 1, N 
D = 0.0E0 

IF (ISW .EQ. 2) GO TO 920 
IF (ALFI(J) .NE. 0.0E0) GO TO 945 

DO 890 I = 1, N 

IF ( ABS ( Z ( I , J) ) .GT. D) D = ABS(Z(I,J)) 
CONTINUE 


C 

C 

_C 

C- 

C 

c 

~c 

c 

c 

_c 

c 

c 

c 

c 

c 

c 

-c 

c 

c 

—C 

c 

w 

c 

- c 

c 

c 

-c 


DO 900 I = 1, N 
900 Z(I,J) = Z(I,J) / D 

GO TO 950 

920 DO 930 I = 1, N 

R = ABS(Z(I, J-l) ) + ABS (Z(I,J) ) 

IF (R .NE. 0.0E0) R = R * SQRT ( ( Z ( I , J-l) /R) **2 
x +(Z(I,J)/R)**2) 

IF (R .GT. D) D = R 
930 CONTINUE 

DO 940 I = 1, N 

Z(I,J-1) = Z(I,J-1) / D 
Z ( I , J) - Z ( I , J) / D 

940 CONTINUE 

945 ISW = 3 - ISW 

950 CONTINUE 

RETURN 
END 

ROUTINE NAME 
FROM EISPACK 


- PF266=RGG 


EISP7491 

EISP7492 

EISP7493 

EISP7494 

EISP7495 

EISP7496 

EISP7497 

EISP7498 

EISP7499 

EISP7500 

EISP7501 

EISP7502 

EISP7503 

EISP7504 

EISP7505 

EISP7506 

EISP7507 

EISP7508 

EISP7509 

EISP7510 

EISP7511 

EISP7512 

EISP7513 

EISP7514 

EISP7515 

EISP7516 

EISP7517 

EISP7518 

EISP7519 

EISP7520 

EISP7521 

RGG 2 


LATEST REVISION 


PURPOSE 


USAGE 

ARGUMENTS NM 


N 


AUGUST 1,1984 

COMPUTER SCIENCES CORP. , HAMPTON, VA. 


- THIS SUBROUTINE CALLS THE RECOMMENDED 
SEQUENCE OF SUBROUTINES FROM THE EIGENSYSTEM 
SUBROUTINE PACKAGE (EISPACK) TO FIND THE 
EIGENVALUES AND EIGENVECTORS (IF DESIRED) FOR 
THE REAL GENERAL GENERALIZED EIGENPROBLEM AX 
= (LAMBDA) BX. 

- CALL RGG(NM,N,A,B,ALFR,ALFI,BETA,MATZ, Z, IERR) 

- ON INPUT NM MUST BE SET TO THE ROW DIMENSION 
OF THE TWO-DIMENSIONAL ARRAY PARAMETERS AS 
DECLARED IN THE CALLING PROGRAM DIMENSION 
STATEMENT . 

- ON INPUT N IS THE ORDER OF THE MATRICES A 
AND B. 

- ON INPUT A CONTAINS A REAL GENERAL MATRIX. 
MUST BE OF DIMENSION NM X N. 


RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 

RGG 


3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 



RGG 

B - ON INPUT B CONTAINS A REAL GENERAL MATRIX. RGG 

MUST BE OF DIMENSION NM X N. RGG 

RGG 

ALFR - ON OUTPUT ALFR CONTAINS THE REAL PART OF THE RGG 
NUMERATORS OF THE EIGENVALUES. RGG 

MUST BE OF DIMENSION N. RGG 

RGG 

ALFI - ON OUTPUT ALFI CONTAINS THE IMAGINARY PART OF RGG 
THE NUMERATORS OF THE EIGENVALUES. RGG 

MUST BE OF DIMENSION N. RGG 

RGG 

BETA - ON OUTPUT BETA CONTAINS THE DENOMINATORS OF RGG 
THE EIGENVALUES, WHICH ARE THUS GIVEN RGG 

BY THE RATIOS (ALFR+I* ALFI) /BETA. RGG 

COMPLEX CONJUGATE PAIRS OF EIGENVALUES APPEAR RGG 
CONSECUTIVELY WITH THE EIGENVALUE RGG 

HAVING THE POSITIVE IMAGINARY PART FIRST. RGG 

MUST BE OF DIMENSION N. RGG 

RGG 

MATZ - ON INPUT MATZ IS AN INTEGER VARIABLE SET EQUALRGG 
TO ZERO IF ONLY EIGENVALUES ARE RGG 

DESIRED. OTHERWISE IT IS SET TO RGG 

ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND RGG 
EIGENVECTORS . RGG 

RGG 

Z - ON OUTPUT Z CONTAINS THE REAL AND IMAGINARY RGG 

PARTS OF THE EIGENVECTORS IF MATZ IS NOT RGG 

ZERO. IF THE J-TH EIGENVALUE IS REAL, THE RGG 
J-TH COLUMN OF Z CONTAINS ITS RGG 

EIGENVECTOR. IF THE J-TH RGG 

EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY RGG 
PART, THE J-TH AND (J+1)-TH RGG 

COLUMNS OF Z CONTAIN THE REAL AND RGG 

IMAGINARY PARTS OF ITS EIGENVECTOR. THE RGG 

CONJUGATE OF THIS VECTOR IS THE RGG 

EIGENVECTOR FOR THE CONJUGATE EIGENVALUE. RGG 

MUST BE OF DIMENSION NM X N. RGG 

RGG 

I ERR - ON OUTPUT I ERR IS AN INTEGER OUTPUT VARIABLE RGG 

SET EQUAL TO AN ERROR COMPLETION CODE RGG 

DESCRIBED IN THE DOCUMENTATION FOR QZIT RGG 

PF261) . THE NORMAL COMPLETION CODE IS ZERO. RGG 

RGG 

RGG 

REQUIRED ROUTINES - PF2 60=QZHES , PF261=QZIT, PF262=QZVAL, PF263=QZVECRGG 

HC3 18=EPSLON RGG 

RGG 

REMARKS 1 . REFERENCES RGG 

RGG 

FROM THE EISPACK PACKAGE OF EIGENSYSTEM ROUTINES. RGG 

RGG 

2. SUBROUTINE RGG IS A DRIVER ROUTINE WHICH CALLS ROUTINESRGG 

QZHES(PF260) , QZIT(PF261), QZVAL(PF262) , AND RGG 

QZVEC(PF263) . RGG 

RGG 


QZHES (PF260) ACCEPTS A PAIR OF REAL GENERAL MATRICESRGG 
AND REDUCES ONE OF THEM TO UPPER HESSENBERG FORM ANDRGG 
THE OTHER TO UPPER TRIANGULAR FORM USING ORTHOGONALRGG 
TRANSFORMATIONS . RGG 


31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 



c 



RGG 

91 

c 


QZIT (PF261) ACCEPTS A PAIR OF REAL MATRICES, 

ONE OFRGG 

92 

c 


THEM IN UPPER HESSENBERG FORM AND THE OTHER IN 

UPPERRGG 

93 

- c 


TRIANGULAR FORM. IT REDUCES THE HESSENBERG MATRIX TORGG 

94 

c 


QUASI -TRIANGULAR FORM USING ORTHOGONAL TRAN S FORMAT I ON SRGG 

95 

c 


WHILE MAINTAINING THE TRIANGULAR FORM OF THE OTHERRGG 

96 

c 


MATRIX. 

RGG 

97 

c 



RGG 

98 

c 


QZVAL (PF2 62 ) ACCEPTS A PAIR OF REAL MATRICES, 

ONE OFRGG 

99 

c 


THEM IN QUASI -TRIANGULAR FORM AND THE OTHER IN 

UPPERRGG 

100 

-*c 


TRIANGULAR FORM. IT REDUCES THE QUAS I -TRI ANGULARRGG 

101 

c 


MATRIX FURTHER, SO THAT ANY REMAINING 2-BY-2 

BLOCKSRGG 

102 

c 


CORRESPOND TO PAIRS OF COMPLEX EIGENVALUES, AND 

RETURN SRGG 

103 



QUANTITIES WHOSE RATIOS GIVE THE GENERALI ZEDRGG 

104 

c 


EIGENVALUES . 

RGG 

105 

c 



RGG 

106 

c 


QZVEC (PF263 ) ACCEPTS A PAIR OF REAL MATRICES, 

ONE OFRGG 

107 

c 


THEM IN QUASI -TRIANGULAR FORM (IN WHICH EACH 

2-BY-2RGG 

108 

c 


BLOCK CORRESPONDS TO A PAIR OF COMPLEX EIGENVALUES) ANDRGG 

109 

c 


THE OTHER IN UPPER TRIANGULAR FORM. IT COMPUTES THERGG 

110 

— c 


EIGENVECTORS OF THE TRIANGULAR PROBLEM AND TRANSFORMSRGG 

111 

c 


THE RESULTS BACK TO THE ORIGINAL COORDINATE SYSTEM. RGG 

112 

c 



RGG 

113 

n 

EXAMPLE : 

RGG 

114 

~c 



RGG 

115 

c 



RGG 

116 

c 

PROGRAM TRGG ( OUTPUT , TAPE6=0UTPUT ) 

RGG 

117 

~c 



RGG 

118 

c 



RGG 

119 

c 

DIMENSION A(5,5),B(5,5) , ALFR ( 5 ) ,ALFI(5) , BETA (5) , Z(5,5) 

RGG 

120 

-C 



RGG 

121 

c 

N = 5 

i 

RGG 

122 

c 

NM = 

5 

RGG 

123 

c 

MATZ 

= 1 

RGG 

124 

c 



RGG 

125 

c 

DATA 

A / 1 0 • ,2. ,3. ,2*1. ,2. ,12. , 1 . ,2. ,1* , 3 . , 1 . ,11. , 

RGG 

126 

c 

* 

l.,-l.,l.,2.,l.,9. ,3*1. , -1. , 1. , 15. / 

RGG 

127 

-c 



RGG 

128 

c 

DATA 

B / 12 . , 1 . , -1 . , 2 . ,2*1. , 14. , 1. ,-l. , 1. ,-l. ,1. , 

RGG 

129 

c 

* 

16. ,-l. ,1. ,2. ,-l. ,-l. ,12. ,-l. ,3*1. ,-l. ,11. / 

RGG 

130 

o 



RGG 

131 

c 

CALL 

RGG (NM, N , A, B , ALFR, ALFI , BETA, MATZ , Z , IERR) 

RGG 

132 




RGG 

133 

n 

WRITE (6 ,99) I ERR 

RGG 

134 

~c 

WRITE (6, 100) ALFR , ALFI , BETA , ( (Z(I, J) ,1=1,5) ,J=1,5) 

RGG 

135 

C99 

FORMAT (1H1,7HIERR = ,14) 

RGG 

136 

C100 

FORMAT ( 1H0 , 7HALFR = /1H ,5(G8.2,2X)/ 

RGG 

137 

~c 

* 

8H0ALFI = /1H , 5 (G8 . 2 , 2X) / 

RGG 

138 

c 

* 

8H0BETA = /1H ,5(G8.2,2X)/ 

RGG 

139 


* 

5H0Z = / 5 ( 1H , 5 (G8 . 2 , 2X) / ) ) 

RGG 

140 

n 

STOP 


RGG 

141 

c 

END 


RGG 

142 

'*'1 



RGG 

143 


OUTPUT : 

RGG 

144 

“c 



RGG 

145 

c 

IERR 

= 0 

RGG 

146 

o 

V-* 

ALFR 

= 

RGG 

147 

~0 

15. 

7.2 16. 10. 8.6 

RGG 

148 

c 

ALFI 

= 

RGG 

149 

n 

0. 

0. 0. 0. 0. 

RGG 

150 



C BETA = RGG 151 

C 9.9 17. 14. 11. 13. RGG 152 

C Z = RGG 153 

-C .26 - . 59E-01 .23 -.30 -1.0 RGG 154 

C -.85 .39 1.0 -.69 .26 RGG 155 

C 1.0 1.0 .85 .88 . 54E-01 RGG 156 

_C -1.0 .83 -.39 .72 -.46 RGG 157 

C -.45 -.84 .65 1.0 -.19E-01 RGG 158 

C RGG 159 

C RGG 160 

SUBROUTINE diverg (NM, N , A, B , ALFR, ALFI , BETA, MATZ , Z , IERR) EISP7 

C 


implicit real*8 (a-h,o-z) 

INTEGER N,NM, IERR, MATZ 

REAL* 8 A (NM, N) , B ( NM , N ) , ALFR(N) , ALFI (N) , BETA (N) ,Z(NM,N) 
LOGICAL TF 


EISP7613 

EISP7614 

EISP7615 

EISP7616 


zero = 0.0e+00 
IF (N .LE. NM) GO TO 10 
IERR = 10 * N 
GO TO 50 

-c 

10 IF (MATZ .NE. 0) GO TO 20 

C FIND EIGENVALUES ONLY 

_ TF = .FALSE. 

CALL QZHES ( NM , N , A , B , TF , Z ) 

CALL QZIT(NM,N,A,B, zero ,TF,Z,IERR) 

CALL QZVAL ( NM , N , A , B , ALFR , ALFI , BETA , TF , Z ) 

GO TO 50 

C FIND BOTH EIGENVALUES AND EIGENVECTORS . . . 

20 TF = .TRUE. 

- CALL QZHES (NM,N, A, B,TF, Z) 

CALL QZIT (NM, N, A, B , zero ,TF,Z,IERR) 

CALL QZVAL (NM,N, A, B, ALFR, ALFI, BETA, TF,Z) 

_ IF (IERR .NE. 0) GO TO 50 

CALL QZVEC(NM,N, A, B, ALFR, ALFI, BETA, Z) 

50 RETURN 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

- END 
subroutine gotoer 
write (6, 10) 

- 10 format ( 'there is an error in calculating subroutine') 

return 


EISP7617 

EISP7618 

EISP7619 

EISP7620 

EISP7621 

EISP7622 

EISP7623 

EISP7624 

EISP7625 

EISP7626 

EISP7627 

EISP7628 

EISP7629 

EISP7630 

EISP7631 

EISP7632 

EISP7633 

EISP7634 

EISP7635 

EISP7636 


end 


EISP7637 


c 



